图像几何失真与局部增强

参考：http://www.bkjia.com/cjjc/971653.html

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);
}