基于CUDA和OpenCV的投影仪梯形校正实现

1. 投影仪梯形校正介绍

        我们在经常使用投影仪时，由于投影仪高度与投影高度不匹配，经常会出现投影仪灯泡网上扬起或者侧着的情形，这种时候投影在幕布或墙上的画面或呈现梯形的形状，对于我们观看投影带来不好的用户体验。现在的大多数厂商都对此加入了数据梯形校正功能，关于梯形校正，您可以参考如下资料：

• 百度百科：梯形校正

• 论文： 《基于图像空间变换和插值运算的投影仪梯形校正法》

2. 基于CUDA和OpenCV的投影仪梯形校正实现

        本文基于CUDA和OpenCV实现了扇形校正，代码中省去了很多校正相关的知识建议读者先补充后再学习代码！

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "stdio.h"
#include <opencv2/highgui.hpp>    // 使用两个工具来读取数据图像
#include <opencv2/imgproc.hpp>
#ifndef __CUDACC__
#define __CUDACC__
#include "cuda_texture_types.h"
#include <iostream>
#include <cstring>
#endif

#define PI 3.14159265358979323846
#define THETA (PI/6)                          // 角度(0,PI]之间
#define R ((img_src.rows-1)/THETA)            // 扇形的内半径

#define DEBUG 1

using namespace std;
using namespace cv;

//texture<unsigned char, cudaTextureType2D, cudaReadModeElementType> textRef;

// 最邻近插值方法
__device__ inline int2 interpolate_neighbor(double x, double y){ return make_int2(roundf(x), roundf(y));}

__global__ void ImgTrans(unsigned char *img_src, unsigned char *img_dst, int rows_dst, 
	                     int cols_dst, int rows_src, int cols_src, double x0, double y0, double r)
{
	unsigned int idx = blockDim.x*blockIdx.x + threadIdx.x;
	unsigned int idy = blockDim.y*blockIdx.y + threadIdx.y;
	unsigned int ii = idy * cols_dst + idx;
	
	if (idx >= cols_dst || idy >= rows_dst)  return;
	//// idx，idy转化为极坐标
	double rou = sqrtf((idx - x0)*(idx - x0) + (idy - y0)*(idy - y0));
	double theta = atanf((idy - y0)/(idx - x0));
		
	// 需要白色填充请注释这一行，不需要白色请放开填充
	if (theta < -THETA / 2.0 || theta > THETA / 2.0  || rou < r || rou > r + cols_src - 1) img_dst[ii] = 255;
	else
	{
		/* 计算新box对应之前的方形坐标像素点 */
		double x = rou - r;
		double y = theta*(rows_src - 1) / THETA + (rows_src - 1) / 2.0;
		//
		int2 pt = interpolate_neighbor(x, y);
		img_dst[ii] = img_src[cols_src*pt.y+pt.x];
	}
}

__global__ void ImgTransBinLinear(unsigned char *img_src, unsigned char *img_dst, int rows_dst,
	                  int cols_dst, int rows_src, int cols_src, double x0, double y0, double r)
{
	unsigned int idx = blockDim.x*blockIdx.x + threadIdx.x;
	unsigned int idy = blockDim.y*blockIdx.y + threadIdx.y;
	unsigned int ii = idy * cols_dst + idx;

	if (idx >= cols_dst || idy >= rows_dst)  return;
	
	//// idx，idy转化为极坐标
	double rou = sqrtf((idx - x0)*(idx - x0) + (idy - y0)*(idy - y0));
	double theta = atanf((idy - y0) / (idx - x0));

	/* 需要白色填充请注释这一行，不需要白色请放开填充 */
	if (theta < -THETA / 2.0 || theta > THETA / 2.0 || rou < r || rou > r + cols_src - 1)
	{
		img_dst[ii] = 255;
	}
	else
	{
		// 计算新box对应之前的方形坐标像素点
		double x = rou - r;
		double y = theta*(rows_src - 1) / THETA + (rows_src - 1) / 2.0;

		double xx = x - floor(x);
		double yy = y - floor(y);   // x,y在每个方格中的局部坐标

		int2 pt = make_int2(floor(x), floor(y));

		img_dst[ii] = img_src[pt.y * cols_src + pt.x] * (1 - xx)*(1 - yy)
			+ img_src[(pt.y + 1)*cols_src + pt.x] * (1 - xx)*yy
			+ img_src[pt.y*cols_src + pt.x] * (1 - yy)*xx
			+ img_src[(pt.y + 1)*cols_src + pt.x + 1] * xx*yy;
	}
}

int main(int argc, char* argv[]){
	Mat img = imread("SrcImages/bModeSlice0001.bin.jpg");
	Mat img_src; cvtColor(img, img_src, COLOR_BGR2GRAY);   // 将图像转化为灰度图,256*330*1

	//Mat img_src(Size(img.cols,img.rows),CV_8UC1, Scalar::all(255));
	imshow("变换前的图形", img_src);

	// 新图像包围的矩形框大小
	int cols = img_src.cols + ceil(R*(1-cos(THETA/2)));
	int rows = (ceil((R+img_src.cols-1)*sin(THETA/2)-0.5)+1) * 2;
	double r = R;

	// 坐标转换对应的(x0,y0),是浮点数，不是整数,且对应变换后的box，而不是变换前
	double y0 = (rows-1)/2.0f, x0 = -R*cos(THETA/2.0);

#if DEBUG
	printf("变换前的图像大小：%dx%d\n", img_src.rows, img_src.cols);
	printf("输入的角度：%f, 半径为:%f \n",THETA*180/PI,R);
	printf("变换后的图像大小: %dx%d\n",rows,cols);
	printf("变换后中心点的坐标: (%0.4f, %0.4f)\n", x0, y0);
#endif

	size_t bts = sizeof(unsigned char)*rows*cols;
	size_t bts_ori = sizeof(unsigned char)*img_src.rows*img_src.cols;
	cudaEvent_t start, stop;
	float elapsed_t = 0.0f;

	//unsigned char *host_src = (unsigned char *)malloc(bts);   // 变换前的图像
	unsigned  char *host_dst = (unsigned  char *)malloc(bts);   // 变换后的图像
	memset(host_dst,0,bts);

	unsigned  char *dev_src; cudaMalloc((void**)&dev_src, bts_ori);
	unsigned  char *dev_dst; cudaMalloc((void**)&dev_dst, bts);

	cudaMemcpy(dev_src, img_src.data, bts_ori, cudaMemcpyHostToDevice);

	cudaEventCreate(&start); cudaEventCreate(&stop);

	//cudaChannelFormatDesc chanelDesc = cudaCreateChannelDesc<unsigned  char>();       // 声明数据类型
	//cudaArray *cuArray;

	//cudaMallocArray(&cuArray, &chanelDesc, img_src.cols, img_src.rows);    // 声明大小为W*H的CUDA数组
	//// 将img_src的数据传输到host_src中去，这里存在坐标变换
	//cudaMemcpyToArray(cuArray, 0, 0, img_src.data, bts, cudaMemcpyHostToDevice);

	//textRef.addressMode[0] = cudaAddressModeClamp;
	//textRef.addressMode[1] = cudaAddressModeClamp;
	//textRef.normalized = false;                      // 是否对纹理坐标归一化
	//textRef.filterMode = cudaFilterModePoint ;      // 纹理的插值模式,还有cudaFilterModeLinear线性插值
	//
	//cudaBindTextureToArray(&textRef,cuArray,&chanelDesc);                  // 绑定纹理，CUDA数组和纹理参考的链接
	
	dim3 block_t(16,16);
	dim3 grid_t(ceil(cols/16.0f),ceil(rows/16.0f));

	cudaEventRecord(start);
	// 最邻近插值法核函数
	//ImgTrans << <grid_t, block_t >> >(dev_src, dev_dst, rows, cols, img_src.rows,img_src.cols, x0, y0,r);
	// 双线性插值核函数
	ImgTransBinLinear << <grid_t, block_t >> >(dev_src, dev_dst, rows, cols, img_src.rows, img_src.cols, x0, y0, r);

	
	cudaEventRecord(stop);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsed_t, start, stop);
	printf("\n======================扇形变换核函数运行时长：%0.4f ms！======================= \n", elapsed_t);

	cudaMemcpy(host_dst, dev_dst, bts, cudaMemcpyDeviceToHost);


	Mat img_dst(rows, cols, CV_8UC1, host_dst);
	imshow("变换后的图形", img_dst);
	waitKey();

	//cudaUnbindTexture(&textRef);	
	//cudaFreeArray(cuArray);

	cudaFree(dev_src);
	cudaFree(dev_dst);
	free(host_dst);

	return 0;
}

3. 运行结果

3.1 运行控制台

3.2 原图

3.3 校正后的扇形图

4. 其它说明

• 本来是想学习纹理内存自带的插值功能和便捷补全功能的，但是试了很多次，纹理内存的使用并没有成功。我也不知道是自己的问题还是说明，最后就自己撸了最邻近和双线性插值的CUDA实现。如果有热心的读者读到这里并且看出我注释代码中使用纹理内存的错误，请留言指正。

• 本文使用OpenCV只是方便图像的显示，所有核心变换和操作函数并没使用OpenCV实现。