图像结合畸变矫正-数字图像处理

/*-------------------------- Author: 叶耀祖 opencv version: opencv343 visual studio 2017 Date:2018/12/18 EMail:yaozuye@zju.edu.cn --------------------------*/ #include <iostream> #include <opencv2/opencv.hpp> #include <cmath> #include <algorithm> # include<string> #define PI 3.14159 #define PHOTO_W 644 #define PHOTO_H 485 #define PHOTO_AREA (644*485) # define PARAM_SIZE 8 using namespace cv; using namespace std; const Point2f CenterPoint(PHOTO_W / 2, PHOTO_H / 2); void int2str(const int &int_temp, string &string_temp) { stringstream stream; stream << int_temp; string_temp = stream.str(); //此处也可以用 stream>>string_temp } //数字图像处理第一次作业源代码 void HK1() { Mat src = Mat(256, 256, CV_8UC3, Scalar(0, 0, 0)); Mat src2 = Mat(640, 640, CV_8UC3, Scalar(0, 0, 0)); int radius = 300; Point center = Point(320, 320); circle(src2, center, radius, Scalar(255,255,255), -1); Mat rectDst; for (int i = 0; i < 256; i++) { for (int j = 0; j < 256; j++) { uchar Hue = uchar(180.0 / 255 * j); uchar Satu = 255-i ; uchar Value = 255; src.at<Vec3b>(i,j) = Vec3b(Hue, Satu, Value); } } for (int i = 0; i < 640; i++) { for (int j = 0; j < 640; j++) { if (src2.at<Vec3b>(i, j) == Vec3b(255, 255, 255)) { int Hue = atan((320 - i)*1.0 / (j - 320)) / PI * 180 + 90; if (j < 320) { Hue += 180; } uchar cvHue = Hue/2; uchar saturation = 1.0*sqrt((320 - i)*(320 - i) + (j - 320)*(j - 320)) / radius * 255; src2.at<Vec3b>(i, j) = Vec3b(cvHue, saturation, 255); } } } cvtColor(src2, src2, COLOR_HSV2RGB); cvtColor(src, src, COLOR_HSV2RGB); imshow("rect", src); imshow("circle", src2); imwrite("rectDst.bmp", src); imwrite("cirleDst.bmp", src2); waitKey(0); } bool absCompare(int a, int b) { return abs(a) < abs(b); } //kernelSize为滤波器核的边长，k为近邻个数 void kNearMean(Mat &src8uc1, Mat &dst, int k = 1, int kernelSize = 3) { int rows = src8uc1.rows; int cols = src8uc1.cols; int size = (kernelSize - 1) / 2; int arrLen = kernelSize * kernelSize - 1; //此处使用int是为了累加和减法不会溢出 int *sortArray = new int[arrLen]; int originValue; for (int i = size; i < rows - size; i++) { for (int j = size; j < cols - size; j++) { int tempIndex = 0; originValue = int(src8uc1.at<uchar>(i, j)); for (int tempi = i - size; tempi <= i + size; tempi++) { for (int tempj = j - size; tempj <= j + size; tempj++) { if (tempi != i || tempj != j) { sortArray[tempIndex] = int(src8uc1.at<uchar>(tempi, tempj)) - int(originValue); tempIndex++; } } } sort(sortArray, sortArray + arrLen, absCompare); int sumTemp = 0; for (int temp = 0; temp < k; temp++) { sumTemp += sortArray[temp]; } int dstValue = originValue + sumTemp*1.0 / k; dst.at<uchar>(i, j) = uchar(dstValue); } } delete[] sortArray; } void meanFilter(Mat &src8uc1, Mat &dst, int kernelSize = 3) { int rows = src8uc1.rows; int cols = src8uc1.cols; int size = (kernelSize - 1) / 2; int arrLen = kernelSize * kernelSize; int originValue; for (int i = size; i < rows - size; i++) { for (int j = size; j < cols - size; j++) { int tempSum = 0; originValue = int(src8uc1.at<uchar>(i, j)); for (int tempi = i - size; tempi <= i + size; tempi++) { for (int tempj = j - size; tempj <= j + size; tempj++) { tempSum += int(src8uc1.at<uchar>(tempi, tempj)); } } tempSum -= originValue; int dstValue = tempSum/(arrLen-1); dst.at<uchar>(i, j) = uchar(dstValue); } } } //PSNR是"peak signal to noise ratio"的缩写，即峰值信噪比 double PSNR(Mat &originSrcUC1,Mat &outputSrcUC1) { int rows = originSrcUC1.rows; int cols = originSrcUC1.cols; double mse = 0; const int MAXI = 255; for (int i = 0; i < rows; i++) { for (int j = 0; j < cols; j++) { mse += pow(double(originSrcUC1.at<uchar>(i, j)) - double(outputSrcUC1.at<uchar>(i, j)), 2); } } mse /= (rows*cols); double psnr = 10.0*log10(pow(MAXI, 2) / mse); return psnr; } //数字图像处理第二次作业源代码 void HK2() { Mat src0 = imread("lena.bmp"); Mat src = imread("lena_noise.bmp"); Mat dst; imshow("src lena", src0); imshow("src lena with noise", src); cout << "--------------------" << endl; cout << "kernel = 3 x 3" << endl; double max_psnr = 0; int max_psnr_k = 1; for (int i = 1; i <= 8; i++) { vector<Mat> BGR; split(src, BGR); kNearMean(BGR[0], BGR[0], i); kNearMean(BGR[1], BGR[1], i); kNearMean(BGR[2], BGR[2], i); merge(BGR, dst); double psnr = PSNR(src0, dst); if (psnr > max_psnr) { max_psnr = psnr; max_psnr_k = i; } cout << "PSNR(k=" << i << "): " << psnr << endl; } vector<Mat> BGR; split(src, BGR); kNearMean(BGR[0], BGR[0], max_psnr_k); kNearMean(BGR[1], BGR[1], max_psnr_k); kNearMean(BGR[2], BGR[2], max_psnr_k); merge(BGR, dst); cout << "max_psnr = " << max_psnr << " k=" << max_psnr_k << endl; imshow("dstImg", dst); waitKey(0); } void pointNormalize(Point2f &srcPoint) { srcPoint -= CenterPoint; srcPoint.x /= (PHOTO_W / 2); srcPoint.y /= (PHOTO_H / 2); } void pointBackNormalize(Point2f &srcPoint) { srcPoint.x *= (PHOTO_W / 2); srcPoint.y *= (PHOTO_H / 2); srcPoint += CenterPoint; } Point2f cordCalc(Point2f oPoint,float param[2][PARAM_SIZE]) { float x = param[0][0] + oPoint.x*param[0][1] + oPoint.y*param[0][2] + oPoint.x*oPoint.y*param[0][3] + oPoint.x*oPoint.x*param[0][4] + oPoint.y*oPoint.y*param[0][5] + oPoint.x*oPoint.y*oPoint.y *param[0][6] + oPoint.x*oPoint.x*oPoint.y *param[0][7]; float y = param[1][0] + oPoint.x*param[1][1] + oPoint.y*param[1][2] + oPoint.x*oPoint.y*param[1][3] + oPoint.x*oPoint.x*param[1][4] + oPoint.y*oPoint.y*param[1][5] + oPoint.x*oPoint.y*oPoint.y *param[1][6] + oPoint.x*oPoint.x*oPoint.y *param[1][7]; return Point2f(x, y); } //loss function 的偏导数，backLoss[0]为偏L偏x'，backLoss[1]为偏L偏y' void lossFuncDiff(Point2f predictPoint, Point2f cPoint, float &backLossX,float &backLossY) { backLossX = -2 * cPoint.x + 2 * predictPoint.x; backLossY = -2 * cPoint.y + 2 * predictPoint.y; } //param_grad为返回的梯度参数 void gradient_calc(vector<Point2f> originPoint, vector<Point2f> currentPoint, vector<Point2f> predicPointVec, float param[2][PARAM_SIZE], float param_grad[2][PARAM_SIZE]) { int vSize = originPoint.size(); for (int i = 0; i < 2; i++) { for (int j = 0; j < PARAM_SIZE; j++) { param_grad[i][j] = 0; } } float backLossX, backLossY; for (int i = 0; i < vSize; i++) { lossFuncDiff(predicPointVec[i], currentPoint[i], backLossX, backLossY); for (int j = 0; j < 2; j++) { if (j == 0) { param_grad[j][0] += backLossX; param_grad[j][1] += backLossX * originPoint[i].x; param_grad[j][2] += backLossX * originPoint[i].y; param_grad[j][3] += backLossX * originPoint[i].x * originPoint[i].y; param_grad[j][4] += backLossX * originPoint[i].x * originPoint[i].x; param_grad[j][5] += backLossX * originPoint[i].y * originPoint[i].y; param_grad[j][6] += backLossX * originPoint[i].y * originPoint[i].y * originPoint[i].x; param_grad[j][7] += backLossX * originPoint[i].x * originPoint[i].x * originPoint[i].y; } else if (j == 1) { param_grad[j][0] += backLossY; param_grad[j][1] += backLossY * originPoint[i].x; param_grad[j][2] += backLossY * originPoint[i].y; param_grad[j][3] += backLossY * originPoint[i].x * originPoint[i].y; param_grad[j][4] += backLossY * originPoint[i].x * originPoint[i].x; param_grad[j][5] += backLossY * originPoint[i].y * originPoint[i].y; param_grad[j][6] += backLossY * originPoint[i].y * originPoint[i].y * originPoint[i].x; param_grad[j][7] += backLossY * originPoint[i].x * originPoint