邊緣偵測

用平面方程去逼近 3x3 影像區塊時，可假設

      f(x,y) = ax + by + c

當中的 (a,b) 表示此區塊的 gradient (可摹想成傾斜向量)
代入 x,y = [-1,0,1] 可得出此區域的代數描述:

     ┌ -a-b+c  -a+c  -a+b+c ┐
     │   -b+c     c     b+c │
     └  a-b+c   a+c   a+b+c ┘

想偵測區塊中一條水平的線，可使用如下的 mask:

     ┌  0   0   0 ┐     ┌  1   1   1 ┐
     │  1   1   1 │ 或  │  0   0   0 │
     └ -1  -1  -1 ┘     └ -1  -1  -1 ┘

後者較佳，因為直行數值為單調增(減)，有助於設計上的最佳化。
再將 mask 改成加權形式，用以增強邊緣點，寫成：

     ┌  β  α  β ┐     ┌  β  0  -β ┐
     │   0   0  0  │ 與  │  α  0  -α │
     └ -β -α -β ┘     └  β  0  -β ┘

將 mask 和區域的代數描述 做 convolution，分別得到一階導數

      2b(α+2β)   與   2a(α+2β)

把他們分別叫做 Gx, Gy，由於此區域之 gradient 量值為

      √(Gx*Gx + Gy*Gy) = 2(α+2β) * √(a*a+b*b)

當 2(α+2β)=1 時，gradient 才符合標準定義 √(a*a+b*b)
於是求解兩根，得 (α,β) = (1/4, 1/8) = ....
將兩根代回 mask 並乘上 8 轉成整數，便得出 Sobel kernel:
  
     ┌  1   2   1 ┐     ┌  1   0  -1 ┐
     │  0   0   0 │ 與  │  2   0  -2 │
     └ -1  -2  -1 ┘     └  1   0  -1 ┘

這程式分別輸出 horizontal 與 vertical mask 運算的結果。

一般對影像邊緣有 zero padding 與 wrapping 與 omit 3 種料理法，
此處採用省略法。無論用的是何種方法，實際上只能處理中央的
(W-2)*(H-2) 區塊。

int sobelH [3*3] = {-1,-2,-1, 0,0,0, 1,2,1};   //平的
int sobelV [3*3] = {-1,0,1, -2,0,2, -1,0,1};   //直的

template <typename T>      
void edge (T* in, T* out, int w, int h, int* m)    //m = kernel
{
    int c;                  //c = color
    out += (w+1);
    in  += (w+1);
    T* end = out+ (w-2)*h-2;            

    while (out < end) {
        c  = m[0]*in[-w-1]; c += m[1]*in[-w]; c += m[2]*in[-w+1];
        c += m[3]*in[  -1]; c += m[4]*in[ 0]; c += m[5]*in[   1];
        c += m[6]*in[ w-1]; c += m[7]*in[ w]; c += m[8]*in[ w+1];
        if (c < 0)   c = -c; //c/=8;
        if (c > 255) c = 255;
        *out++ = c;
        in++;
    }
}