Geometry & Recognition

기준 좌표계에 대해서 원점을 이동하고 좌표축을 회전시킨 새로운 좌표계에서 점의 좌표는 바뀐다. 원래의 좌표와 바뀐 좌표값 사이의 관계를 주는 변환이 Isometric transformation (isometry)이다. 평면에서 이 변환은 평행이동을 나타내는 파라미터 2개, 그리고 1개의 회전각 파라미터에 의해서 결정이 된다. 회전각이 $θ$고, 평행이동이 $(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) \rightarrow(u_1, v_1), (x_2, y_2)\rightarrow(u_2, v_2)\}$ 만 있으면 이들 파라미터를 정확히 결정할 수 있다. 그러나 변환에 필요한 점 정보를 얻는 과정은 필연적으로 노이즈의 영향을 받게 되므로 주어진 모든 점을 정확히 연결하는 변환을 일반적으로 구할 수 없다. 이 경우에는 isometry 파라미터는 일반적으로 최소자승법에 의해서 결정될 수 있다. 

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

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

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

$$L = \sum_i \big [ (a  x_i - b  y_i - u_i)^2 + (b  x_i + a  y_i - v_i)^2 \big] + λ  (a^2 + b^2 - 1) ;$$

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

$$\sum_i  ( a  x_i - b  y_i - u_i) x_i + ( b  x_i + a  y_i - v_i) y_i + λ a = 0 $$

$$ \sum_i  ( a  x_i - b  y_i - u_i) (-y_i) + ( b  x_i + a  y_i - v_i) x_i + λ b = 0$$

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

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

$$a = ∑(x_i u_ i + y_ i v_ i) / ∑ (x_ i^2 + y_i^2 + λ)$$

$$ b = ∑ (x_i v_ i - y_i u_i) / ∑ (x_i^2 + y_i^2 + λ) $$
또한, Lagrange 멀티플라이어 $λ$는

$$A  = ∑ (x_i u_i + y_i v_i)$$

$$B =  ∑ (x_i v_i - y_i u_i);$$

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

$$∑ ( x_i^2 + y_i^2 + λ ) = \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) \rightarrow (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개와 회전각($θ$) 1개, 그리고 $x, y$축으로의 평행이동을 나타내는 parameter ($t_x$, $t_y$) 2 개를 합해서 총 4개가 있어야 한다. 이 parameter에 의해서 원본 이미지의 픽셀 $(x, y)$가 변환된 이미지의 픽셀 $(u, v)$에 대응한다고 하면, 이들 간의 관계는 다음식으로 주어진다.

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

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

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

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

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

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

\begin{align}L = \sum_{i} &\Big[ \left| u_i - \left(s\cos(θ) x_i - s \sin(θ) y_i + t_x\right)\right|^2 \\ &+ \left|v_i - \left(s \sin(θ) x_i + s \cos(θ) y_i + t_y\right)\right|^2\Big]\end{align}

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

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

$$ \frac {\partial L}{\partial b}=0:\quad \sum _{i} (u_i - (ax_i - by_i + t_x))(y_i) + (v_i - (bx_i + a y_i + t_y))(-x_i) = 0$$

$$\frac {\partial L}{\partial t_x}=0: \quad \sum_{i} (u_i - (ax_i - by_i + t_x)) = 0$$

$$ \frac {\partial L}{\partial t_y}=0: \quad  \sum_{i} (vi - (bx_i + ay_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_iy_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_iy_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 + bS_y - t_x  N = 0$$

$$ S_v - b S_x - aS_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})}\begin{bmatrix}  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{bmatrix} $$

아래의 코드는 이것을 구현한 것이다. 물론, $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) \rightarrow (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

이미지에서 관찰된 점집합이 $\{(x_i, y_i)| i = 1, 2,\dots, N\}$이 있다. 이 점집합을 직선 $y = a + bx$ 로 피팅을 하고 싶을 때, 보통 최소자승법을 이용하는데, 원리는 직선의 방정식이 예측한 $y$값과 실제 관찰한 $y$값의 차이의 제곱(=square deviation)을 최소화시키는 직선의 기울기 $a$와 절편 $b$를 찾는 것이다:

$$\chi^2(a, b) = \sum_i |y_i - (b x_i +a) |^2 $$

데이터를 얻는 측정 과정에서 측정값 $y_i$는 랜덤 노이즈를 포함하게 되고, 이 노이즈는 참값 $y(x)$ 근방에서 정규분포를 한다고 가정을 할 수 있다. 만약 모든 측정의 노이즈가 동일한 표준편차 $\sigma$를 갖게 된다면, $N$개의 관측 데이터가 나타날 확률(likelihood)은 (개별 측정은 IID 조건을 만족한다고 전제)

$$P = \prod_{i} e^{ -\frac{ | y_i -  (bx_i + a)|^2 }{2\sigma^2 }  }$$

의 형태가 된다. 따라서 최소자승법은 이 likelihood를 최대화시키는 파라미터를 찾는 방법이 된다. 최소자승법은 피팅 파라미터를 주어진 데이터를 이용해서 표현할 수 있는 장점은 있지만, outliers에 대해서 매우 취약하다 (아래의 결과 그림을 참조). 이는 적은 수의 outlier도 χ2에 큰 기여를 할 수 있기 때문이다. 따라서 피팅을 좀 더 robust 하게 만들기 위해서는 outliers가 likelihood에 기여하는 정도를 줄여야 한다. 이를 위해서는 likelihood의 지수 인자를 큰 에러에서 덜 민감하게 반응하는 꼴로 바뀌어야 한다. 이를 만족하는 가장 간단한 것 방법 중 하나가 square-deviation 대신에 absolute-deviation을 이용하는 것이다:

$$\text{absolute deviation} = \sum _i | y_i - (bx_i + a)|   .$$

그러나 이 식을 사용하는 경우에는 최소자승법과 다르게 기울기 $a$와 절편 $b$를 주어진 데이터 $\{(x_i, y_i)\}$로 표현할 수 없고, 반복적인 방법을 이용해서 찾아야 한다. 

수열 $\{c_i\}$에 대해 합 $\sum_{i} |c_i - a|$은 $a$가 수열의 median 값일 때 최솟값을 갖는다는 사실을 이용하면 (증명: 극값을 구하기 위해서 $a$에 대해서 미분하면, $0=(\sum_{c_i > a} 1)-(\sum_{c_i < a} 1)$: 합은 $a$가 $c_i$ 보다 큰 경우와 작은 경우로 분리. 따라서 0이 되기 위해서는 작은 경우와 큰 경우의 수가 같아야 한다. 고로, $a = \text{median}\{c_i\}$ q.e.d.). 고정된 절편 $b$에 대해서 absolute deviation을 최소로 만드는 기울기 $a$는

$$a= \text{median} \{ y_i - b x_i\}$$

임을 알 수 있다. 그리고  absolute deviation 식을 절편 $b$에 대해서 미분해서

$$0 = \sum_i \text{sign} \left( y_i - (bx_i +a) \right)$$

을 얻는데, 위에서 구한 기울기 $a$를 대입한 후 bracketing and bisection 방법을 이용하면 절편 $b$를 얻을 수 있다(불연속 함수이므로 일반적으로 근을 구하는 방법을 사용하는 것은 위험하다). 아래의 코드는 이를 구현한 것이다.

double FitLine_LS(std::vector<double>& x, std::vector<double>& y, double *a, double *b);

// 최소자승법을 이용한 직선 추정:
// return (sigma[dy] / sigma[x]);
double FitLine_LS(std::vector<double>& x, std::vector<double>& y, double& a, double& b) {
    double sx = 0, sy = 0, sxx = 0, sxy = 0;
    for (int i = x.size(); i-->0;) {
        sx  += x[i];        sy  += y[i];
        sxx += x[i] * x[i]; sxy += x[i] * y[i];
    };
    const int n = x.size();
    double det = n * sxx - sx * sx;
    if (det == 0.) return -1;                   // vertical line;
    a = (sxx * sy - sx * sxy) / det;
    b = (n * sxy - sx * sy) / det;
    double chi2 = 0;
    for (int i = x.size(); i-->0;) {
        double t = y[i] - (*a + *b * x[i]);
        chi2 += t * t;
    }
    det /= n;         //det -> var(x) * n;
    // chi2 = var(dy) * n;
    // (dy vs x의 편차비)
    return  sqrt(chi2 / det);
}

// 기울기(bb)가 주어진 경우에 y-절편(median = aa)값을 추정하고, 이 때 AD값을 얻는다.
double RhoFunc(std::vector<double>& x, std::vector<double>& y,
               double bb, double& aa, double& abdev) {
    std::vector<double> h(x.size());
    for (int i = x.size(); i-->0;) h[i] = y[i] - bb * x[i];
    std::sort(h.begin(), h.end());
    // median;
    const int med = h.size()/2;
    aa = (h.size() & 1) ? h[med] : (h[med] + h[med-1])/2;

    double sum = 0;
    abdev = 0;
    for (int i = x.size(); i-->0;) {
        double d = y[i] - (bb * x[i] + aa);
        abdev += fabs(d);
        // sign-함수의 원점에서 모호함을 없애기 위해서 증폭을 시킴;
        if (y[i] != 0.) d /= fabs(y[i]);
        if (fabs(d) > DBL_EPSILON) // sign 함수의 모호함에서 벗어나면,
            sum += (d >= 0 ? x[i] : -x[i]);
    }
    return sum; // return sum{xi * sign(yi - (b * xi + a))}
};
// y = a + b * x ;
// Least Absolute Deviation:
double FitLine_MAD (std::vector<double>& x, std::vector<double>& y,
                    double& a, double& b) {
    // least square esimates for (aa, bb);
    double aa, bb, abdev;
    double sigb = FitLine_LS(x, y, aa, bb);   // estimate: y=aa + bb*x;
    double b1 = bb;
    double f1 = RhoFunc(x, y, b1, aa, abdev);
    /* bracket 3-sigma away in the downhill direction;*/
    double b2 = fabs(3 * sigb);
    b2 = bb + (f1 < 0 ? -b2 : b2);
    double f2 = RhoFunc(x, y, b2, aa, abdev);

    // if conditional added to take care of the case of a
    // line input into this function. It is needed to avoid an
    // infinite loop when (b1 == b2) (within floating point error)
    if (fabs(b2 - b1) > (sigb + 0.005)) {
        // bracketing;
        while ((f1 * f2) > 0) {
            bb = 2 * b2 - b1;
            b1 = b2; b2 = bb; 
            f1 = f2; f2 = RhoFunc(x, y, b2, aa, abdev) ;
        }
    }
    // refine until the error is a negligible number of std;
    sigb *= 0.01;
    while (fabs(b2 - b1)> sigb) {
        // bisection;
        bb = (b1 + b2) / 2.;
        if ((bb == b1) || (bb == b2)) break ;
        double f = RhoFunc(x, y, bb, aa, abdev) ;
        if ((f * f1) >= 0) {
            f1 = f; b1 = bb;
        } else {
            f2 = f; b2 = bb;
        }
    }
    a = aa; b = bb; 
    return (abdev/x.size());
}

// 붉은 선--> 최소자승법에 의한 피팅.: outlier에 매우 취약한 구조.
// 파란 선--> least absolute deviation을 이용한 피팅: outlier에 매우 robust 하다.