数字图象处理——仿射变换

在《数字图象处理》中提供了基于

图像仿射矩阵T。图像旋转、偏移变换等变换原理相同，根据这些在自己实现时，可分为三个步骤：

1. 坐标转换
   这些变换是以图像的中心为坐标原点（O’(x’,y’)），而图片原本设定是以图像左上角为坐标原点（O(x,y)）。所以要进行坐标的变换。
   
   用矩阵的方式表示：
   
   逆运算：
   

2. 图像变换
   变换公式：
   
   逆运算：
   

3. 坐标还原
   因为此时，坐标原点是在图像的中心，为方便操作需要将坐标原点重新变换到图像左上角。
   从旋转后到旋转前的坐标变换为：
   
   （w’，h’是旋转后的图像的宽高）；
   而逆运算为：
   （1）

下面是旋转的最邻近内插法和双线性内插法的实现：
其实在实现中，最重要的是求出变换后的图像宽高，再根据（1）式推出变换后的每一个像素点对应的原图像坐标。根据（1）式计算原图像像素点坐标x，y时，结果中分别会有一个常量，用dx，dy表示。

void nearestInterpolation(Mat &src, Mat &dst, float dx, float dy, double theta);
int mainqqw()
{
    cv::Mat src = imread("D:/xitong/picture/rain.jpg");
    namedWindow("orginal");
    imshow("orginal", src);
    int srcwidth = src.cols;
    int srcheigh = src.rows;
    /*旋转角度*/
    double theta = 30.0f*3.1415926 / 180.0f;
    /*
        转换坐标原点到图像中心
                ∧y
                |
            0   |   1
                |
        --------o--------->x
                |
            2   |   3
                |
    */
    float srcX[4], srcY[4];
    srcX[0] = (float)(-((srcwidth - 1) / 2));
    srcX[1] = (float)((srcwidth - 1) / 2);
    srcX[2] = (float)(-(srcwidth - 1) / 2);
    srcX[3] = (float)((srcwidth - 1) / 2);
    srcY[0] = (float)((srcheigh - 1) / 2);
    srcY[1] = (float)((srcheigh - 1) / 2);
    srcY[2] = (float)(-(srcheigh - 1) / 2);
    srcY[3] = (float)(-(srcheigh - 1) / 2);

    /*
        旋转后的图像坐标，此时坐标原点依然是旋转中心
    */
    float dstX[4], dstY[4];
    for (int i = 0; i < 4; i++)
    {
        dstX[i] = cos(theta)*srcX[i] + sin(theta)*srcY[i];
        dstY[i] = -sin(theta)*srcX[i] + cos(theta)*srcY[i];
    }

    /*
        ==>旋转后图像长宽
    */
    int dstwidth = (max(fabs(dstX[3] - dstX[0]), fabs(dstX[2] - dstX[1])) + 0.5);
    int dstheigh = (max(fabs(dstY[3] - dstY[0]), fabs(dstY[2] - dstY[1])) + 0.5);

    /*Mat dst = Mat(Size(src.rows * 2, src.cols * 2), src.type(), Scalar::all(0));
    nearestInterpolation(src, dst, 0.5);*/
    Mat dst;
    dst.create(dstheigh, dstwidth, src.type());

    //（1）式在运算后的常量，为运算方便提前的出结果
    float dx = -0.5*dstwidth*cos(theta) - 0.5*dstheigh*sin(theta) + 0.5*srcwidth;
    float dy = 0.5*dstwidth*sin(theta) - 0.5*dstheigh*cos(theta) + 0.5*srcheigh;

    nearestInterpolation(src, dst, dx, dy, theta);
    waitKey();
    return 0;
}

/*
    最邻近内插旋转
*/
void nearestInterpolation(Mat &src, Mat &dst, float dx, float dy, double theta)
{
    int x, y;
    for (int i = 0; i < dst.rows; i++)
    {
        for (int j = 0; j < dst.cols; j++)
        {
            /*
                根据（1）式，推出对应原图像 像素坐标的运算结果
            */
            x = cvFloor(float(j)*cos(theta) + float(i)*sin(theta) + dx);
            y = cvFloor(float(-j)*sin(theta) + float(i)*cos(theta) + dy);
            if ((x < 0) || (x >= src.cols) || (y < 0) || (y >= src.rows))
            {
                if (src.channels() == 3)
                    dst.at<Vec3b>(i, j) = Vec3b(0, 0, 0);
                if (src.channels() == 1)
                    dst.at<uchar>(i, j) = 0;
            }
            else
            {
                if (src.channels() == 3)
                    dst.at<Vec3b>(i, j) = src.at<Vec3b>(y, x);
                if (src.channels() == 1)
                    dst.at<uchar>(i, j) = src.at<uchar>(y, x);
            }

        }
    }
    namedWindow("最邻近内插旋转");
    imshow("最邻近内插旋转", dst);
}
/*
    双线性内插法旋转
*/
void bilinearRotate(Mat &src, Mat &dst, float dx, float dy, double theta)
{
    float fu, fv;
    int x, y;
    Vec3b point[4];
    uchar upoint[4];
    for (int j = 0; j < dst.rows; j++)
    {
        for (int i = 0; i < dst.cols; i++)
        {
            fu = float(j)*cos(theta) + float(i)*sin(theta) + dx;
            fv = float(-j)*sin(theta) + float(i)*cos(theta) + dy;
            x = cvFloor(fu);
            y = cvFloor(fv);
            fu -= x;
            fv -= y;

            if ((x < 0) || (x >= src.cols-1) || (y < 0) || (y >= src.rows-1))
            {
                if (src.channels() == 3)
                    dst.at<Vec3b>(i, j) = Vec3b(0, 0, 0);
                if (src.channels() == 1)
                    dst.at<uchar>(i, j) = 0;
            }
            else
            {
                if (src.channels() == 3)
                {
                    point[0] = src.at<Vec3b>(y, x);
                    point[1] = src.at<Vec3b>(y + 1, x);
                    point[2] = src.at<Vec3b>(y, x + 1);
                    point[3] = src.at<Vec3b>(y + 1, x + 1);
                    dst.at<Vec3b>(i, j) = (1 - fu)*(1 - fv)*point[0] + (1 - fu)*(fv)*point[1] + (1 - fv)*(fu)*point[2] + fu*fv*point[3];
                }

                if (src.channels() == 1)
                {
                    upoint[0] = src.at<uchar>(y, x);
                    upoint[1] = src.at<uchar>(y + 1, x);
                    upoint[2] = src.at<uchar>(y, x + 1);
                    upoint[3] = src.at<uchar>(y + 1, x + 1);
                    dst.at<uchar>(i, j) = (1 - fu)*(1 - fv)*upoint[0] + (1 - fu)*(fv)*upoint[1] + (1 - fv)*(fu)*upoint[2] + fu*fv*upoint[3];
                }
            }
        }
    }
    namedWindow("双线性内插法旋转");
    imshow("双线性内插法旋转", dst);
}