Geometry & Recognition

Harris response:

$$ R=\text{Det}(M)-k\text{Trace}(M)^{2}~\quad k \in [0.04,0.06]$$

참고: http://www.ipol.im/pub/art/2018/229/article_lr.pdf

double PI = 4.0 * atan(1.0);
double Gaussian(double x, double y, double sigma) {
    double sigma2 = sigma * sigma;
    double t = (x * x + y * y) / (2 * sigma2);
    double n = 1.0 / (2.0 * PI * sigma2);
    return n * exp( -t );
};

struct Corner {
    double x, y;
    double strength;
    int    flag;
};
int HarrisCornerDetector(BYTE *img, int w, int h,
                    double sigma/*=1.2*/, double kvalue/*=0.04~0.06*/,
                    double minMeasure/*=1e-1*/,
                    std::vector<Corner>& corners)
{
    // Sobel gradient 3x3
    const int nn[] = { -w - 1, -w, -w + 1, -1, 0, 1, w - 1, w, w + 1}; // 3x3 window
    const int sobelX[] = { -1, 0, 1, -2, 0, 2, -1, 0, 1};
    const int sobelY[] = { -1, -2, -1, 0, 0, 0, 1, 2, 1};
    const int xmax = w - 1, ymax = h - 1;
    std::vector<double> Gx(w * h, 0);
    std::vector<double> Gy(w * h, 0);
    for (int y = 1, pos = w; y < ymax; y++) { // pos = w; starting address;
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            double sx = 0, sy = 0;
            for (int k = 0; k < 9; k++) {
                int v = img[pos + nn[k]];
                sx += sobelX[k] * v;
                sy += sobelY[k] * v;
            }
            Gx[pos] = sx / (4 * 255);
            Gy[pos] = sy / (4 * 255);
        }
        pos++; // skip x=xmax;
    }
    // precompute the coefficients of the Gaussian filter
    const int radius = int(2 * sigma);
    const int window = 1 + 2 * radius;
    std::vector<double> gFilter(window * window);
    for (int j = -radius, pos = 0; j <= radius; j++)
        for (int i = -radius; i <= radius; i++, pos++)
            gFilter[pos] = Gaussian(i, j, sigma);

    // compute the Harris measure for each pixel of the image
    std::vector<double> harrismap(w * h, 0);
    double hmax = 0;
    for (int y = 0, pos = 0; y < h; y++) {
        for (int x = 0; x < w; x++, pos++) {
            // Convolve gradient with Gaussian filter:
            double filtersum = 0, sxx = 0, syy = 0, sxy = 0;
            for (int yk = y - radius, wpos = 0; yk <= y + radius; yk++) {
                for (int xk = x - radius; xk <= x + radius; xk++, wpos++) {
                    if (yk < 0 || yk >= h) continue;
                    if (xk < 0 || xk >= w) continue;
                    // Gaussian weight
                    double gauss = gFilter[wpos];
                    filtersum += gauss;
                    // convolution;
                    int loc = yk * w + xk; //image location;
                    sxx += gauss * Gx[loc] * Gx[loc];
                    syy += gauss * Gy[loc] * Gy[loc];
                    sxy += gauss * Gx[loc] * Gy[loc];
                }
            }
            sxx /= filtersum;
            syy /= filtersum;
            sxy /= filtersum;
            // Harris corner measure = Det(M) - kvalue.(Trace(M))^2
            double hmeasure = sxx * syy - sxy * sxy - kvalue * (sxx + syy) * (sxx + syy);
            if (hmeasure <= 0) continue;
            if (hmeasure > hmax) hmax = hmeasure;
            harrismap[pos] = hmeasure;
        }
    }
    // log scale
    double scale = 255 / log(1 + hmax);
    for (int pos = w * h; pos-- > 0;)
        harrismap[pos] = scale * log(1 + harrismap[pos]);
    // Nonmaximal suppression: keep only local maxima;
    minMeasure *= 255; //255 = the maximum value of Harris measure(in log scale);
    const int nn8[] = { -w - 1, -w, -w + 1, 1, w + 1, w, w - 1, -1}; // 8-neighbors;
    for (int y = 1, pos = w; y < ymax; y++) {
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            // thresholding: Harris measure > epsilon
            double hmeasure = harrismap[pos];
            if (hmeasure <= minMeasure) continue;
            bool isLocalMax = true;
            for (int k = 0; k < 8; k++) {
                double hk = harrismap[pos + nn8[k]];
                if (hk >= hmeasure) {
                    isLocalMax = false;
                    break;
                }
            }
            if (isLocalMax) {
                Corner c;
                c.x = x, c.y = y, c.strength = hmeasure, c.flag = 1;
                corners.push_back(c);
            }
        }
        pos++; //skip x=w-1;
    }

    // remove corners too close to each other: keep the highest measure
    double minDistance = 8;
    std::vector<Corner>::iterator iter = corners.begin();
    while (iter != corners.end()) {
        for (int i = 0; i < corners.size(); i++) {
            Corner& c1 = corners[i];
            if (c1.flag == 0) continue;
            if ((c1.x == iter->x && c1.y == iter->y) ) continue;
            if (sqrt((c1.x-iter->x)*(c1.x-iter->x) + (c1.y-iter->y)*(c1.y-iter->y)) > minDistance) 
                continue;
            if (c1.strength < iter->strength) continue;
            iter->flag = 0;    // iter = corners.erase(iter); berased = true;
            break;
        }
        iter++;
    }
    return 1;
};

Otsu 이진화에서 threshold 값의 선택은 전경과 배경 사이의 분산을 최대로 하는 픽셀 값이다. 영상의 히스토그램이 bimodal인 경우는 잘 동작하지만, unimodal인 영상의 경우는 제대로 처리가 안된다. Valley emphasis 방법은 Otsu 방법을 사용하면서 histogram profile의 valley에 가중치를 더해주는 방식으로 threshold 값을 찾는다.

\begin{gather} \text{Ostu's method: } ~~~\text{threshold}= \underset{0\le t< 256}{\text{argmax}}~\left(\omega_1(t) m_1(t)^2 + \omega_2(t) m_2^2(t)\right)\\  \text{valley-emphasis: }~~\text{threshold}=\underset{0\le t < 256} {\text{argmax}}~  (1-p_t)\left(\omega_1(t) m_1(t)^2 + \omega_2(t) m_2^2(t)\right) \end{gather}

Bimodal 분포의 영상은 Otsu와 거의 같은 결과를 준다. 두 개 이상의 클래스로 분리도 쉽게 구현이 된다.

참고: http://msn.iecs.fcu.edu.tw/report/data/ori_paper/2006-12-21/2006-Automatic%20thresholding%20for%20defect%20detection.pdf

int valley_emphasis_Otsu(int hist[], int N) {
    int istart = 0, iend = N - 1;
    while (!hist[istart]) istart++;
    while (!hist[iend])   iend--;
    double st = 0, sxt = 0;           // total pdf, total cdf;
    for (int i = istart; i <= iend;i++) { 
        st += hist[i]; 
        sxt += double(i) * hist[i];
    }
    int winner = istart; 
    double maxgain = 0;
    double s0 = 0, sx0 = 0;
    for(int i = istart; i <= iend; i++) {
        s0  += hist[i];				// 1-pdf;
        sx0 += double(i) * hist[i];		// 1-cdf;
        double s1  = st - s0;			// 2-pdf
        double sx1 = sxt - sx0;			// 2-cdf
        double m0  = sx0 / s0;			// E(X1);
        double m1  = sx1 / s1;			// E(X2);
        double gain = (1 - hist[i] / st) * (s0 * m0 * m0 + s1 * m1 * m1);
        if (gain > maxgain) {
            maxgain = gain; 
            winner = i;
        }
    }
    return winner;
};

평면을 한 점 $(x, y)$이 원점을 축으로 하는 회전에 의해서 $(x', y')$ 위치로 옮겼다면 두 점 사이의 관계는 

$$ \begin{pmatrix} x' \\ y'\end{pmatrix} = \begin{pmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta\end{pmatrix}\begin{pmatrix} x\\ y\end{pmatrix} = \mathbf{R}_\theta \begin{pmatrix} x\\ y\end{pmatrix}$$

로 표현된다. 여기서 $\theta$는 시계방향으로 회전한 각이다. 이 회전 변환은 3번의 shear 변환의 곱으로 다음처럼 쓸 수 있다.

$$ \mathbf{R}_\theta = \mathbf{S}_x. \mathbf{S}_y . \mathbf{S}_x$$

$$ \mathbf{S}_x =\begin{pmatrix} 1 & \tan (\theta/2) \\0 &1\end{pmatrix}, ~~\mathbf{S}_y =\begin{pmatrix}1 & 0 \\ -\sin \theta & 1\end{pmatrix}$$

이제 Fourier변환과 shear 변환 사이의 관계를 알아보자. 일반적인 2차원 Fourier 변환은

\begin{gather} F(u, v)  =\iint   f(x, y) e^{-2\pi i (ux + vy)}dxdy={\cal F}[ f(x, y)] \\= {\cal F_x}\left[{\cal F_y} [f(x, y)]\right]={\cal F_y}\left[  {\cal  F_x}  [f(x, y)]\right] = F_1 (u) F_2 (v)\end{gather}

여기서 $\cal F_{x}, F_y$는 각각 $x$, $y$-방향 Fourier 변환이고 순서에 상관이 없다.

이제 $f(x, y)$에 $x$-방향으로 $a$ 만큼 shear 변환을 한 결과가 $g(x, y)$일 때,

$$ g(x, y) = f(x + ay, y)$$

로 쓸 수 있고, 여기에 $x$ 방향으로 Fourier 변환을 하면

$$ G_1(u, y) = {\cal F_x}[g(x, y)]= {\cal F_x}[f(x+ay, y)] = F_1(u, y) e^{-2\pi i ua y} =e^{-2\pi i  ua y}  {\cal F_x}[ f(x, y)]$$ 

임을 알 수 있다. 즉, shear 변환된 함수의 Fourier 변환은 원함수의 Fourier변환에 shear 파라미터에 의존하는 위상 요소 $e^{-2\pi i  ua y}$만 곱해주면 쉽게 구할 수 있음을 알 수 있다. 또한 shear 변환에 의해서 스펙트럼이 위상만 변화고 크기에는 변화가 생기지 않음도 알 수 있다.

참고: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.190.4473&rep=rep1&type=pdf

아래의 코드는 이미지의 중앙을 기준으로 회전을 한 경우를 보여준다: $\theta = 60^\circ$. 중간 단계에서 artifact를 없애기 위해 원본 이미지의 가장자리를 $0$으로 채워 두 배 크기로 만든 후 수행한다. 임의의 지점에 대한 회전도 shear 변환 코드를 조금 손보면 쉽게 구현할 수 있다.

중간단계: 

int fft1d(int dir, int nn, double x[], double y[]); /* https://kipl.tistory.com/22 */
#define HORIZONTAL 0
#define VERTICAL   1
/* perform shear transform along "axis" direction; */
void shearImage(int axis, double shear,
                int width, int height, double *input, double *output) 
{
    int size, halfsz, othersize, halfothersz;
    double TWOPI = 8.0 * atan(1.);
    /* size is the size of the image in the direction of the axis of shear. 
    ** othersize the size of the image in the orthogonal direction 
    */
    if (axis == HORIZONTAL) {
        halfsz = (size = width) >> 1;
        halfothersz = (othersize = height) >> 1;
    }
    else {
        halfsz = (size = height) >> 1;
        halfothersz = (othersize = width) >> 1;
    }

    std::vector<double> re_sig(size, 0);
    std::vector<double> im_sig(size, 0);
    for (int i = 0; i < othersize; i++){
        if (axis == HORIZONTAL)
            memcpy(&re_sig[0],  &input[i * size], size * sizeof(double));
        else
            for (int j = 0; j < size; j++) re_sig[j] = output[j * othersize + i];

        std::fill(im_sig.begin(), im_sig.end(), 0);  // im: padding 0;
        fft1d(1, size, &re_sig[0], &im_sig[0]);
        double shift = - (i - halfothersz - .5) * shear  * TWOPI;
        // No need to touch j = 0;
        for (int j = 1; j < halfsz; j++) {
            double cc = cos(j * shift / size);
            double ss = sin(j * shift / size);
            double reval = re_sig[j];
            double imval = im_sig[j];
            re_sig[j] = cc * reval - ss * imval;
            im_sig[j] = ss * reval + cc * imval;
            re_sig[size - j] =  re_sig[j];
            im_sig[size - j] = -im_sig[j];
        }
        re_sig[halfsz] = cos(shift/2) * re_sig[halfsz] - sin(shift/2) * im_sig[halfsz];
        im_sig[halfsz] = 0.;

        fft1d(-1, size, &re_sig[0], &im_sig[0]);
        if (axis == HORIZONTAL)
            memcpy(&output[i * size], &re_sig[0], size * sizeof(double));
        else
            for (int j = 0; j < size; j++) output[j * othersize + i] = re_sig[j];
    }
};

주어진 함수 $y(x)$를 구간 $0\le x < L$에서 일정한 간격으로 샘플링해서 얻은 $N$개의 data $\{y_n = y(nL/N)|n=0, 1,...,N-1\} $의 discrete Fourier transform(DFT) $Y_k$는

$$ Y_ k = \frac {1}{N} \sum _{n = 0}^{N-1} y_n e^{ -\frac {2\pi i}{N} n k}, $$

로 정의된다. 역변환(IDFT)은 $Y_k$에서 $y_n$을 얻는 과정으로 다음 식으로 계산된다:

$$ y_n = \sum_{k=0}^{N-1} Y_k e^{+\frac {2\pi i}{N}  k n}. $$

이산적인 data $y_n$만을 가지고 $x_n = nL/N$에서 미분계수를 구하려면 어떻게 해야 할까? 샘플링 점 사이 구간에서 값을 추정할 수 있도록 보간 함수를 만들어 사용해야 한다. 간단하게는 IDFT 식에서 $n$번째 샘플링 위치인 $x_n = n(L/N)$을 실수 $x$로 대체하여 연속함수를 얻는 방법을 생각할 수 있다. 그리고 IDFT 식에서 $ Y_k e^{+\frac {2\pi i}{N} k n}$ 대신에 추가적으로 위상 요소를 곱해서 $Y_k e^{+\frac {2\pi i}{N} (k + m_k N) n}, ~m_k =\text {integer}$ 로 바꾸더라도 동일한 원 데이터 $y_n$을 얻는다는 사실을 이용하도록 하자. 이 두 가지 관찰을 이용해서 만든 보간 함수는

$$ y(x) = \sum _{k=0}^{N-1} Y_k e^{+\frac{2\pi i}{L} (k +m_k N ) x}  $$

이다. 추가 위상이 $y_n$은 바꾸지 않지만 어떤 선택을 하는지에 따라 샘플링 사이 구간에서 다른 진동을 보이므로 미분계수가 영향을 받게 된다. 그러면 어떤 정수 $m_k$을 선택해야 하는가? 여러 가지 선택이 있을 수 있지만 전구간에서 도함수의 제곱 평균을 최소로 만드는 $m_k$을 선택하도록 하자: 

$$\underset {m_k} {\text {argmin} } ~ \frac {1}{L} \int_0^L | y'(x)|^2 dx  $$

그런데 $$\frac{1}{L} \int_0^L | y'(x)|^2 dx  = \Big(\frac {2\pi}{L} \Big)^2 ~\sum_{k=0}^{N-1}  (k+m_k N)^2|Y_k|^2  $$ 이므로 각 $k$에 대해서 $|k +m_k N|^2$이 최소가 되도록 선택해야 한다. $k<N/2$일 때는 $m_k=0$을 선택하면 되고, $k>N/2$일 때는 $m_k=-1$을 선택하면 된다. $N=\text{even}$일 때는 $m_{N/2}=0 \text{ or }1$이 같은 결과를 주는데, 두 기여를 균등하게 되도록 선택한다: 

$$ Y_{N/2} e^{\frac{2\pi i}{L} (N/2 + m_{N/2} N)x} \rightarrow \frac{1}{2} Y_{N/2}\left( e^{\frac{\pi i}{L} Nx} + e^{-\frac{\pi i}{L} Nx} \right) = Y_{N/2} \cos \Big( \frac{\pi}{L} Nx\Big) $$ 따라서 $m_k$ 선택은 정리하면 다음과 같다:

$$  m_k= \left\{\begin{array}{ll}   0, & 0\le k < N/2   \\ -1, &  N/2 <k <N  \end{array} \right.$$

따라서 주어진 샘플링 데이터 $\{y_n\}$의 DFT 계수 $\{ Y_n\}$을 구하면 샘플링 데이터를 보간하는 함수(Fourier interpolation)는

$$ y(x) = Y_0 + \sum_{0 <k <N/2}  \left( Y_k e^{+\frac {2\pi i}{L}  k x} + Y_{N-k} e^{ -\frac {2\pi i}{L}  k x} \right)  +  Y_{N/2} \cos\Big( \frac{\pi}{L}Nx\Big) $$

로 얻어진다(마지막 항은 $N$이 짝수일 때만 더해짐). 이 보간함수의 도함수는 

$$y'(x)= \sum_{0 <k <N/2} \frac {2\pi i }{L} k \Big(Y_k e^{+\frac {2\pi i}{L}  k x} - Y_{N-k} e^{ -\frac{2\pi i}{L}  k x}  \Big) - \frac{\pi}{L} N Y_{N/2} \sin \Big( \frac{\pi}{L} Nx\Big)$$

로 계산된다(마지막 항은 $N$이 짝수일 때만 더해짐). 샘플링 위치에서 미분계수는

$$ y'_n = y'(nL/N) = \sum_{0<k<N/2} \frac{2\pi i}{L} k \Big( Y_k e^{+\frac{2\pi i}{N} kn} - Y_{N-k} e^{ -\frac{2\pi i}{N} kn}  \Big) = \sum_{k=0}^{N-1} Y'_k e^{+\frac{2\pi i}{N} kn  }$$

처럼 쓸 수 있는데 이 식은 $Y'_k$의 IDFT 형태이다. 

따라서 DFT를 이용해서 미분계수를 구하는 전체적인 알고리즘은

1. DFT로 $Y_k$을 구하고,

2. $Y'_k$는 $$  Y_k' = \left\{  \begin {array}{ll}    \frac {2\pi i}{L} k  \times Y_k & 0 \le k <N/2\\ \frac{2\pi i}{L} (k-N) \times Y_k &  k > N/2 \\ 0 & k=N/2 ~\text {when } N=\text {even} \end {array} \right.  $$

3. $Y'_k$에 IDFT를 해서 $y'_n$을 구한다.

위에서 구한 보간 함수를 이용하면 이차 미분 값도 DFT를 이용해서 구할 수 있다.

아래 그림은 6개의 주파수 $f_k = 2,5,9,11,21,29\text{Hz}$가 섞인 신호 함수의 미분을 FFT를 이용해서 구한 결과를 정확한 도함수(green)와 함께 보여준다.

$$f(x) = \sum_{k=0}^{5} (k+1) \cos (2\pi f_k x)$$

$0\le x < 1=L$ 구간에서 일정한 간격으로 샘플링한 $N=64$개의 데이터를 사용하였다. 

void fftgrad1d() {
    const double TWOPI = 8.0 * atan(1.0);
    const int samples = 64;
    double re[samples], im[samples];
    const int nfreqs = 6;
    const double freq[nfreqs] = {2, 5, 9, 11, 21, 29};
    for ( int i = 0; i < samples; i++ ) {
        double ti = double(i) / samples;
        double signal = 0;
        for ( int k = 0; k < nfreqs; k++ )
            signal += (k + 1) * cos ( TWOPI * freq[k] * ti);
        re[i] = signal;
        im[i] = 0.;
    }
    // int fft1d(int dir, int nn, double x[], double y[]);
    fft1d(1, samples, re, im);
    const double norm = TWOPI / 1.0;
    for (int i = 0; i < samples; i++) {
        double rv = re[i];
        double iv = im[i];
        double fac = norm * double (i > samples / 2 ? i - samples: i);  
        if (i == samples / 2) fac = 0;
        re[i] = - fac * iv;
        im[i] = fac * rv;
    }
    //invers FFT: re[] stores derivatives.
    fft1d(-1, samples, re, im);
}

이차원에서 DFT를 이용해서 gradient을 구하는 경우는 복소함수 $f = Re(x, y) + i Im(x, y)$의 gradient가

$$\Big(\frac{\partial}{ \partial x} +i \frac{\partial}{\partial y} \Big) (Re + i Im)= i\Big( \frac{\partial Re}{\partial y} + \frac{\partial Im}{\partial x} \Big) - \Big(- \frac{\partial Re}{\partial x} + \frac{\partial Im}{\partial y}\Big)  $$

임을 이용하면 된다.

영상(실수 함수)의 IDFT 식에 적용하면 미분 연산자가 우변에 전체적으로 $i$를 곱하는 효과를 주게 되고, 우변식의 실수부와 허수가 각각 영상의 $x$-방향/ $y$-방향 미분에 해당한다.

void FFTGrad(CRaster& raster, CRaster& out) {
    const double TWOPI = 8.0 * atan(1.0);
    CSize sz = raster.GetSize();
    int nx = sz.cx, ny = sz.cy;
    int stride = nx << 1;
    std::vector<double> data(stride * sz.cy);
    for (int y = 0; y < ny; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        double *re = &data[y * stride];
        double *im = re + nx;
        for (int x = 0; x < nx; x++) {
            re[x] = *p++;
            im[x] = 0;
        }
    }
    // forward FFT(1): https://kipl.tistory.com/22
    fft2d(&data[0], nx, ny, 1);
    double normx = TWOPI / double(nx);
    double normy = TWOPI / double(ny);
    int hnx = nx >> 1;
    int hny = ny >> 1;
    for (int y = 0; y < ny; y++) {
        double *re = &data[y * stride];
        double *im = re + nx;
        for (int x = 0; x < nx; x++) {
            double rv = re[x];
            double iv = im[x];
            double xx = normx * double(x > hnx ? x - nx: x);
            double yy = normy * double(y > hny ? y - ny: y);
            if (x == hnx) xx = 0;  // N = even;
            if (y == hny) yy = 0;
            re[x] =  yy * rv + xx * iv;
            im[x] = -xx * rv + yy * iv;
        }
    }
    /* inverse FFT of {Y'_k} */
    fft2d(&data[0], nx, ny, -1);
    std::vector<double> mag(nx * ny);
    double magmax = 0, magmin = 1.e+20;
    for (int y = 0, k = 0; y < ny; y++) {
        double *re = &data[y * stride];
        double *im = re + nx;
        for (int x = 0; x < nx; x++) {
            double m = hypot(re[x], im[x]);
            if (m > magmax) magmax = m;
            if (m < magmin) magmin = m;
            mag[k++] = m;
        }
    }
    // stretch the magnitude;
    for (int y = 0, k = 0; y < ny; y++)
        for (int x = 0; x < nx; x++) {
            double a = 255 * (mag[k++] - magmin) / (magmax - magmin);
            out.SetPixel(x, y, 255-int(a));
        }
}

최소 자승법 문제는 추정치와 실제의 측정치와의 차이를 최소로 만드는 parameter vector를 구하는 것이다.
$$\underset{\mathbf{x}} {\text{argmin}} ~|\mathbf{A}. \mathbf{x} - \mathbf{b}|^2.$$

여기서 design matrix $\mathbf{A}$는 추정에 사용이 된 basis 함수를 각각의 독립변수에서 계산한 값이고, $\mathbf{x}$는 구하고자 하는 parameter vector이며, $\mathbf{b}$는 측정값 vector이다. 예를 들면 주어진 측정값을 $n-1$차의 다항식을 이용하여서 피팅하려고 하는 경우에는 parameter는 다항식의 계수가 되며 $n$-차원의 vector space을 형성하게 되며, $\mathbf{A}$는

$$   A_{ij} = (x_i)^j ,  \quad j=0,..., n-1$$

가 일 것이다. 일반적으로 $A_{ij}$는 $x_i$에서 계산된 basis-함수의 값이 된다. 위의 식을 $\mathbf{x}$에 대해서 미분을 하면 극값을 취하면 최소자승 문제는 아래의 행렬식을 푸는 문제가 된다

$$    (\mathbf{A}^T. \mathbf{A}) .\mathbf{x} =  \mathbf{A}^T.  \mathbf{b}.$$

$\mathbf{A}^T. \mathbf{A}$ 은 $n\times n$ matrix다. 이 행렬이 역행렬을 가지게 되면

 $$ \mathbf{x} = (\mathbf{A}^T. \mathbf{A})^{-1} . (\mathbf{A}. \mathbf{b}),$$

를 하여서 원하는 parameter vector를 얻을 수 있다. 그러나 피팅 문제에 따라 행렬 $\mathbf{A}^T. \mathbf{A}$가 매우 singular 해져 역행렬을 구할 수 없게 되는 경우에 종종 생긴다. 예를 들면, 저주파의 신호를 고주파 기저 함수를 이용하여서 최소자승법을 사용하고자 하는 경우 등에 이러한 문제에 부딪히게 된다. 이런 경우에는 직접적으로 $\mathbf{A}^T. \mathbf{A}$의 역행렬을 구하는 방법을 이용할 수 없고

$$   \mathbf{A} .\mathbf{x} =  \mathbf{b}$$

의 식을 $\mathbf{A}$의 SVD(Singular Value Decomposition)를 이용하여서 풀 수가 있다. $\mathbf{A}$를 SVD 하면 $\mathbf{A}_{m\times n}=\mathbf{U}_{m\times n} . \mathbf{w}_{n\times n}. \mathbf{V}_{n\times n}^T $의 형태로 분해할 수 있다. 여기서 $\mathbf{w}=\text{diag}(\underbrace{w_0, w_1,...}_{\text{nonzero}},0,..,0)$로 쓰여지는 대각행렬이다. matrix $\mathbf{U}$와 $\mathbf{V}$의 column vector를 사용하면
$$ \mathbf{A}  =\sum_{w_k \ne 0} w_k \mathbf{u}_k \otimes \mathbf{v}_k^T$$

의 형태로 쓰인다. $\mathbf{u}_k$는 $\mathbf{U}$의 $k$-번째 열벡터이고, $\mathbf{v}_k$는 $\mathbf{V}$의 $k$-번째 열벡터로 각각 orthonormal basis를 형성한다. parameter 벡터를 $\{ \mathbf{v}_k \}$ basis로 전개를 하면 영이 아닌 singularvalue에 해당하는 성분만 가지게 된다. 구체적으로 위의 $\mathbf{A}$ 분해와 $\mathbf{u}_j^T.\mathbf{u}_k=\delta_{jk}$, 그리고 $\sum_k \mathbf{v}_k \otimes \mathbf{v}_k^T= \mathbf{I}_{n\times n}$임을 이용하면,

\begin{gather}  \mathbf{v}_k^T . \mathbf{x} = \mathbf{u}_k^T . \mathbf{b} / w_k, \quad w_k \ne 0, \\                    \mathbf{v}_k^T . \mathbf{x} = 0, \quad w_k = 0, \\  \rightarrow ~~\mathbf{x} = \sum _{w_k \ne 0 } ( \mathbf{u}_k^T . \mathbf{b} / w_k)  \mathbf{ v} _k , \end{gather}

이어서 위의 해를 구할 수 있다. 이 해는 $|\mathbf{A} . \mathbf{x} -  \mathbf{b}|^2$를 최소화한다.

int svd(double *A, int m, int n, double* w, double *V); // from cmath libary.
void fit_func(double x, double val[], int n) {          // polynomial fitting sample;
    val[0] = 1;
    for(int i = 1; i < n; ++i)
        val[i] = x * val[i - 1];
}
#define EPSILON 1.E-8
int svd_fit(const double x[], const double y[], const int m, const int n,
            void (*fit_func)(double , double [], int ),
            double params[],
            double *error)
{
    double *A = new double [m * n];
    double *w = new double [n];
    double *V = new double [n * n];
    // evaluate design matrix;
    for (int i = 0; i < m; ++i)
        fit_func(x[i], &A[i * n + 0], n) ;

    svd(A, m, n, w, V);
    // now A becomes U matrix;
    // truncate small singular values;
    double wmax = 0;
    for (int i = 0; i < n; ++i)
        if (w[i] > wmax) wmax = w[i];
    double thresh = wmax * EPSILON;
    for (int i = 0; i < n; ++i)
        if (w[i] < thresh) w[i] = 0;
    
    // back substitution;
    double *tmp = new double [n];
    for (int j = 0; j < n; ++j) {
        double s = 0;
        if (w[j]) {
            for (int i = 0; i < m; ++i)
                s += A[i * n + j] * y[i];
            s /= w[j];
        }
        tmp[j] = s;
    }
    for (int j = 0; j < n; ++j) {
        double s = 0;
        for (int jj = 0; jj < n; ++jj)
            s += V[j * n + jj] * tmp[jj];
        params[j] = s;
    };

    //estimate error;
    *error = 0;
    for (int i = 0; i < m; ++i) {
        fit_func(x[i], &A[i * n + 0], n); //use A as a tmp buffer;
        double sum = 0;
        for (int j = 0; j < n; ++j) sum += params[j] * A[i * n + j] ;
        double err = (y[i] - sum);
        *error += err * err ;
    }
    delete[] A; delete[] w; delete[] V;
    delete[] tmp;
    return 1;
}