【OpenCV】幾何失真校正與局部增強，opencv幾何失真校正

【OpenCV】幾何失真校正與局部增強，opencv幾何失真校正

1. 原理

圖像在獲取過程中，由於成像系統的非線性、飛行器姿態的變化等原因，成像後的圖像與原景物圖像相比，會產生比例失調，甚至扭曲。這類圖像退化現象稱之為幾何失真（畸變）。
產生這種原因有：成像系統本身具有的非線性，攝像時視角的變化，被攝對象表面彎曲等。例如，由於視像管攝像機及陰極射線管顯示器的掃描偏轉系統有一定的非線性，常常枕形失真或者桶形失真；由於斜視角度獲得的圖像透視失真等等。
幾何失真主要是由於圖像中的像素點發生位移而產生的，其典型表現為圖像中的物體扭曲、遠近比例不協調等。解決這類失真問題的方法成為幾何畸變校正，簡稱為幾何校正。
有成像系統引起的幾何失真的校正方法有兩種：一種是預畸變法，即采用與畸變相反的非線性掃描偏轉法，用來抵消預計的圖像畸變；另一種方法是所謂的後驗校正法，使用多項式曲線在水平和垂直方向去擬合每一畸變的網線，然後求的反變換的校正函數，用這個校正函數即可校正畸變圖像。
幾何畸變校正分為兩步：第一步是對原圖像坐標空間進行幾何變換，以使像素落在正確的位置上；第二步是重新確定新像素的灰度值，這是因為經過上面的坐標變換後，有些像素點有時被擠壓在一起，有時又被分散開，使校正後的像素不落在離散的坐標點上，因此需要重新確定這些像素的灰度值。 

• 幾何畸變的描述

任意幾何級那都可以由非失真坐標系（x,y）變換到失真坐標系（x',y'）的方程來定義。
設f(x,y)是無失真的原始圖像，g(x',y')是f(x,y)畸變的結果，這一失真的過程是已知的，並且可用函數h1(x,y)和h2(x,y)定義，於是有:
這是幾何校正的基本關系式，這種失真的復原問題實際上是映射變換問題。 

• 幾何變換

從幾何校正的基本關系可見，已知畸變圖像g(x',y')的情況下要求原始圖像f(x,y)的關鍵是要求的函數h1(x,y)和h2(x,y)，則f(x,y)的求取方法就較為簡單了。但實際中往往h1(x,y)和h2(x,y)不知道，這時我們可以采用後驗校正法。
通常h1（x,y）、h2(x,y)可用多項式來近似式中，N為多項式的次數，aij、bij為各自項的待定系數。
後驗校正方法的思想是通過一些已知的正確像素點和畸變點之間的對應關系，擬合出以上兩個多項式的系數，擬合出的多項式作為恢復其他畸變點的變換基礎。例如，一個基准圖通過成像系統後形成畸變圖像，通過研究基准圖像與畸變圖像之間的對應關系，找出多項式的各系數。
N=1時，變換是線性的通常也可以用這種線性畸變來近似較小的幾何畸變。然而由於實際情況復雜多樣，上式在N=2以上就不一定有解惑找不到最優解了，這時就要用最小二乘法了。 

• 局部增強

局部增強方法是根據所關心的局部區域的特性來計算變換或濾波參數，並將結果用於所關心的局部區域，以得到所需的相應的增強效果。由此可見，局部增強方法比全局增強方法在具體進行增強操作前多了一個選擇確定局部區域的步驟，而對每個局部區域仍可采用前幾節介紹的增強方法進行增強。
直方圖變換是空域增強中最常采用的方法，它也很容易用於圖像的局部增強。只需先將圖像分成一系列（一般互相重疊的）小區域（子圖像），此時直方圖均衡化或規定化都可以基於小區域內的像素分布進行，從而使各小區域得到不同的增強效果。

2. 方法

采用內插法確定像素的灰度值。當原圖像坐標（x,y）變換後，落在畸變圖像內，但不是剛好在圖像像素點上，就需要通過一定的手段求出這點的灰度值，常用的方法有最近鄰法，雙線性內插法和三次卷積法。 

• 最近鄰法

較簡單的插值方法是最近鄰法，及選擇離它所映射到的位置最近的輸入像素的灰度值為差值結果。若原圖像上坐標為（x,y）的像素經變換後落在畸變圖像g(其中這種插值法對於鄰近像素點的灰度值有較大改變，但細微結構是粗糙的。 

• 雙線性內插法

原圖像f(x,y)上的一像素坐標為(x,y)，經變換後，落在畸變圖像g(x',y')內的坐標為(u,v)，下面式中[]表示取整。
定義：a = u-[u],b = v-[v],則g(u,v)的取值按如下公式計算
當u=[u]或v=[v]時，則有
復原圖像這就是雙線性內插法。
與最近鄰法相比，內插法幾何校正灰度連續，結果一般滿足要求，但計算量較大且具有低通特性，圖像輪廓模糊。如果要進一步改善圖像質量，可以選用三次卷積法。 

• 局部增強

局部增強也可在對整幅圖增強時直接利用局部信息以達到不同局部不同增強的目的。例如局部空域增強中有一種常用的方法是利用每個像素的鄰域內像素的均值和方差這兩個特性進行的。這裡均值是一個平均亮度的測度，而方差是一個反差的測度。具體來說，如要把輸入圖f(x,y)增強成輸出圖g(x,y) ，需要在每個像素位置(x,y)進行如下變換：
稱為局部增益函數。
以上兩式中，m（x,y)是以像素（x,y)為中心的鄰域內的灰度均值；σ（x,y)是以像素（x,y)為中心的鄰域內的灰度均方差值； M是以f（x,y) 的平均灰度值；k 是一個比例常數。
A（x,y)與 f（x,y)和 m（x,y)的差相乘能放大圖像的局部變化，因為 A（x,y)反比於均方差，所以在圖像中對比度較小的區域得到的增益反而較大，這樣就可以取得局部增強的效果。
式（4.7.1）中最後又將 m（x,y)加回去是為了恢復原區域的平均灰度值。實際中為了平衡圖像中孤立區域灰度值的偏移，常常只將 m（x,y)的一部分加回去，並常將 A（x,y）限制在一定的范圍內。

3. 關鍵代碼

• 最近鄰插值（前向映射）

int main()
{
    IplImage* img = cvLoadImage("Image/slant.jpg",0);
    IplImage* out = cvCreateImage(cvGetSize(img), img->depth, img->nChannels);

    //手動指定四邊形頂點,A、B、C、D為映射前（傾斜），Ao、Bo、Co、Do為映射後（校正）
    CvPoint A = cvPoint(145,1);
    CvPoint B = cvPoint(475,191);
    CvPoint C = cvPoint(355,399);
    CvPoint D = cvPoint(25,209);

    CvPoint Ao = cvPoint(60,80);
    CvPoint Bo = cvPoint(440,80);
    CvPoint Co = cvPoint(440,320);
    CvPoint Do = cvPoint(60,320);

    float m[16] = {    A.x, A.y, A.x*A.y, 1,
                    B.x, B.y, B.x*B.y, 1,
                    C.x, C.y, C.x*C.y, 1,
                    D.x, D.y, D.x*D.y, 1};
    
    float n1[4] = { Ao.x, Bo.x,    Co.x, Do.x};
    float n2[4] = { Ao.y, Bo.y,    Co.y, Do.y};

    CvMat* M = cvCreateMat(4,4,CV_32FC1);
    CvMat* N1 = cvCreateMat(4,1,CV_32FC1);
    CvMat* N2 = cvCreateMat(4,1,CV_32FC1);
    CvMat* K4 = cvCreateMat(4,1,CV_32FC1);
    CvMat* K8 = cvCreateMat(4,1,CV_32FC1);

    cvSetData(M, m, CV_AUTOSTEP);
    cvSetData(N1, n1, CV_AUTOSTEP);
    cvSetData(N2, n2, CV_AUTOSTEP);

    //求解前向映射方程組,取得8個參數
    cvSolve(M, N1, K4, CV_LU);
    cvSolve(M, N2, K8, CV_LU);

    float k1 = K4->data.fl[0];
    float k2 = K4->data.fl[1];
    float k3 = K4->data.fl[2];
    float k4 = K4->data.fl[3];
    float k5 = K8->data.fl[0];
    float k6 = K8->data.fl[1];
    float k7 = K8->data.fl[2];
    float k8 = K8->data.fl[3];

    CvScalar s;    //source

    //前向映射
    for(int x=0; x<img->width; x++)
    {
        for(int y=0; y<img->height; y++)
        {            
            //最近鄰
            int i = int(k1*x + k2*y + k3*x*y + k4 +0.5);
            int j = int(k5*x + k6*y + k7*x*y + k8 +0.5);

            if(i<0 || i>=out->width || j<0 || j>= out->height)
            {
                //旋轉後超出邊緣部分丟棄，空白部分默認用白色填補
            }
            else
            {
                s = cvGet2D(img,y,x);
                cvSet2D(out,j,i,s);
            }
        }
    }
    
    cvNamedWindow("Source", CV_WINDOW_AUTOSIZE);
    cvNamedWindow("Front", CV_WINDOW_AUTOSIZE);
    cvShowImage("Source", img);
    cvShowImage("Front", out);
    cvWaitKey(0);
    cvReleaseImage(&img);
    cvReleaseImage(&out);
    cvDestroyWindow("Source");
    cvDestroyWindow("Front");
} 

• 最近鄰插值（後向映射）

與前向映射的區別之處：

    //手動指定四邊形頂點,A、B、C、D為映射前（校正），Ao、Bo、Co、Do為映射後（傾斜）
    CvPoint A = cvPoint(60,80);
    CvPoint B = cvPoint(440,80);
    CvPoint C = cvPoint(440,320);
    CvPoint D = cvPoint(60,320);

    CvPoint Ao = cvPoint(145,1);
    CvPoint Bo = cvPoint(475,191);
    CvPoint Co = cvPoint(355,399);
    CvPoint Do = cvPoint(25,209);

……
    
    //後向映射
    for(int x=0; x<out->width; x++)
    {
        for(int y=0; y<out->height; y++)
        {    
            //最近鄰
            int i = int(k1*x + k2*y + k3*x*y + k4 +0.5);
            int j = int(k5*x + k6*y + k7*x*y + k8 +0.5);

            if(i<0 || i>=img->width || j<0 || j>= img->height)
            {
                //旋轉前在邊緣外部分默認用白色填補
            }
            else
            {
                s = cvGet2D(img,j,i);
                cvSet2D(out,y,x,s);
            }
        }
    }

…… 

• 雙線性插值（後向映射）

與最近鄰插值區別之處：

……
CvScalar s,a,b,c,d;  //a,b,c,d 為所求s的4個最近鄰像素

    //後向映射
    for(int x=0; x<out->width; x++)
    {
        for(int y=0; y<out->height; y++)
        {    
            //雙線性
            float xx = k1*x + k2*y + k3*x*y + k4;
            float yy = k5*x + k6*y + k7*x*y + k8;
            int i = int(xx);
            int j = int(yy);

            if(i<0 || i>=img->width-1 || j<0 || j>= img->height-1)
            {
                //旋轉前在邊緣外部分默認用白色填補,邊界上的像素這裡不處理
            }
            else
            {
                a = cvGet2D(img,j ,i);
                b = cvGet2D(img,j ,i+1);
                c = cvGet2D(img,j+1 ,i);
                d = cvGet2D(img,j+1 ,i+1);

                float e = (xx - i)*(b.val[0] - a.val[0]) + a.val[0];
                float f = (xx - i)*(d.val[0] - c.val[0]) + c.val[0];

//計算插入的灰度值
                int value = int((yy - j)*(f - e) + e + 0.5);
                s = cvScalar(value,0,0);
                cvSet2D(out,y,x,s);
            }
        }
}
…… 

• 局部增強

IplImage* img = cvLoadImage("Image/leaf.jpg",0);
    IplImage* out = cvCreateImage(cvGetSize(img), img->depth, img->nChannels);    
    float k = 0.5;    //增強系數
    
    CvScalar avg = cvAvg(img);
    float M = avg.val[0]; //計算平均灰度值M    
    std::pair<float, float> cal;
    for(int i=1;i<img->width-1;i++)
{
        for(int j=1;j<img->height-1;j++)
        {
            cal = Calculate(img,i,j);  //自定義方法
            int cur = cvGet2D(img,j,i).val[0];
            int enh;
            
            if( cal.second == 0 )    //防止除數為0;如不處理，在計算增強後的灰度時會導致結果為0，圖像出現黑斑。

                enh =cur;
            else
                //利用局部增強函數計算灰度值並四捨五入
                enh = k * (M / cal.second) * (cur - cal.first) + cal.first + 0.5;
            
            //灰度值裁剪,避免無效灰度值
            enh = enh<0 ? 0 : enh;
            enh = enh>255 ? 255 : enh;
            
            cvSet2D(out,j,i,cvScalar(enh,0,0));
        }
}
……
//自定義的方法,求3*3鄰域的均值與標准差
std::pair<float,float> Calculate(IplImage* img, int i, int j)
{
    int sum = 0;
    int data[9];    
    int n = 0;
    for(int a=-1;a<2;a++)
    {
        for(int b=-1;b<2;b++)
        {
            data[n] = cvGet2D(img,j+b,i+a).val[0];
            sum += data[n];
            n++;
        }
}
    float avg = sum/9.0;
    float tmp = 0;
    
    for(n=0;n<9;n++)
        tmp += pow(data[n]-avg,2);
    
    float sd = sqrt(tmp/9.0);
    return std::make_pair(avg,sd);
}

4. 結果

配置好OpenCV與Visual Studio2012環境後執行上述代碼，輸出結果如下圖所示。(為了排版，縮小了圖片大小。請適當放大觀察，效果更好。)
Figure 1 原始圖片
 
Figure 2 最近鄰插值 
（前向映射，出現空隙並且顆粒明顯）
 
Figure 3 最近鄰插值 
（後向映射，無空隙但顆粒依然明顯）
 
Figure 4 雙線性插值 
（後向映射，無空隙並且過渡也相對自然）

Figure 5 原圖
（葉脈較為模糊，細節不分明）
 
Figure 6 局部增強後
（增強系數k=0.5，紋理更為清晰）