基于C++&CUDA实现全景图转换为天空盒图像

基于C++在CPU中运行：

#include     
#include           
#include    
#include
#include 

#define OutSize 1000
using namespace cv;
using namespace std;
float M_PI = 3.14159265358979323846f;


float faceTransform[6][2] =
{
	{ 0, 0 },  //front
	{ M_PI / 2,0 },//right
	{ M_PI,0 }, //back
	{ -M_PI / 2,0 },//left
	{ 0,-M_PI / 2 },//top
	{ 0,M_PI / 2 }//bottom
};


inline void createCubeMapFace(const Mat &in, Mat &face, int faceId = 0, const int width = -1, const int height = -1)
{
	float inWidth = in.cols;
	float inHeight = in.rows;     // 获取图片的行列数量
	Mat mapx(height, width, CV_32F);
	Mat mapy(height, width, CV_32F);                                         //分配图的x,y轴
																			 // Calculate adjacent (ak) and opposite (an) of the
																			 // triangle that is spanned from the sphere center 
																			 //to our cube face.
	const float an = sin(M_PI / 4);
	const float ak = cos(M_PI / 4);                                          //计算相邻ak和相反an的三角形张成球体中心												
	const float ftu = faceTransform[faceId][0];
	const float ftv = faceTransform[faceId][1];
	// For each point in the target image, 
	// calculate the corresponding source coordinates.                      对于每个图像计算相应的源坐标
	for (int y = 0; y < height; y++) {
		for (int x = 0; x < width; x++) {
			// Map face pixel coordinates to [-1, 1] on plane               将坐标映射在平面上
			float nx = (float)y / (float)height - 0.5f;
			float ny = (float)x / (float)width - 0.5f;
			nx *= 2;
			ny *= 2;
			// Map [-1, 1] plane coord to [-an, an]                          
			// thats the coordinates in respect to a unit sphere 
			// that contains our box. 
			nx *= an;
			ny *= an;
			float u, v;
			// Project from plane to sphere surface.
			if (ftv == 0) {
				// Center faces
				u = atan2(nx, ak);
				v = atan2(ny * cos(u), ak);
				u += ftu;
			}
			else if (ftv > 0) {
				// Bottom face 
				float d = sqrt(nx * nx + ny * ny);
				v = M_PI / 2 - atan2(d, ak);
				u = atan2(ny, nx);
			}
			else {
				// Top face
				//cout << "aaa";
				float d = sqrt(nx * nx + ny * ny);
				v = -M_PI / 2 + atan2(d, ak);
				u = atan2(-ny, nx);
			}
			// Map from angular coordinates to [-1, 1], respectively.
			u = u / (M_PI);
			v = v / (M_PI / 2);
			// Warp around, if our coordinates are out of bounds. 
			while (v < -1) {
				v += 2;
				u += 1;
			}
			while (v > 1) {
				v -= 2;
				u += 1;
			}

			while (u < -1) {
				u += 2;
			}
			while (u > 1) {
				u -= 2;
			}

			// Map from [-1, 1] to in texture space
			u = u / 2.0f + 0.5f;
			v = v / 2.0f + 0.5f;

			u = u*(inWidth - 1);
			v = v*(inHeight - 1);

			mapx.at(x, y) = u;
			mapy.at(x, y) = v;
		}
	}

	// Recreate output image if it has wrong size or type. 
	if (face.cols != width || face.rows != height ||
		face.type() != in.type()) {
		face = Mat(width, height, in.type());
	}

	// Do actual  using OpenCV's remap
	remap(in, face, mapx, mapy, CV_INTER_LINEAR, BORDER_CONSTANT, Scalar(0, 0, 0));
	if (faceId == 0)
	{
		imwrite("cube0000.jpg", face);
	}
	else if (faceId == 1)
	{
		imwrite("cube0001.jpg", face);
	}
	else if (faceId == 2)
	{
		imwrite("cube0002.jpg", face);
	}
	else if (faceId == 3)
	{
		imwrite("cube0003.jpg", face);
	}
	else if (faceId == 4)
	{
		imwrite("cube0004.jpg", face);
	}
	else if (faceId == 5)
	{
		imwrite("cube0005.jpg", face);
	}
	

}


void main() {
	int width = OutSize;
	int height = OutSize;
	cv::Mat srcimage = cv::imread("input.jpg");
	cv::Mat resultImage;
	for (int i = 0; i < 6; i++) {
		createCubeMapFace(srcimage, resultImage, i, width, height);
	}
}

基于CUDA的GPU运算：

createCubeMapFace(gpu_panorama, cubeImages, diff_area, Index_x, Index_y, cubeout);

	//上述参数定义
 	cv::cuda::GpuMat gpu_warp;
	gpu_warp.upload(host_srcimg);
	cv::cuda::GpuMat cubeImages = cv::cuda::GpuMat(cube_W, cube_W, device_srcimg.type());
	cv::cuda::GpuMat diff_area = cv::cuda::GpuMat(pano_H, pano_W, CV_32FC1);
	cv::cuda::GpuMat Index_x = cv::cuda::GpuMat(pano_H, pano_W, CV_32FC1);
	cv::cuda::GpuMat Index_y = cv::cuda::GpuMat(pano_H, pano_W, CV_32FC1);
	//diff_area是索引全景图坐标映射到哪个天空面；
	//Index_x, Index_y是为了定位在天空面的坐标索引。
	//上面两个索引功能不需要的话自动忽略。

__global__ void kernel7(cv::cuda::PtrStepSz de_img, HOD_Point2f dev_faceTransform[], int* i, float* dev_an, float* dev_ak, cv::cuda::PtrStepSzf mapx, cv::cuda::PtrStepSzf mapy, cv::cuda::PtrStepSzf diff_area, cv::cuda::PtrStepSzf Index_x, cv::cuda::PtrStepSzf Index_y)
{
	const int X = blockIdx.x * blockDim.x + threadIdx.x;
	const int Y = blockIdx.y * blockDim.y + threadIdx.y;


	float ftu = dev_faceTransform[(*i)].x;
	float ftv = dev_faceTransform[(*i)].y;
	float u, v;

	if (X < mapy.cols && Y < mapy.rows)
	{

		float nx = ((float)Y / (float)mapx.rows - 0.5) * 2.0 * (*dev_an);
		float ny = ((float)X / (float)mapx.cols - 0.5) * 2.0 * (*dev_an);

		if (ftv == 0) {
			// Center faces
			u = dev_atan2f(nx, (*dev_ak));
			v = dev_atan2f(ny * dev_cos(u), (*dev_ak));
			u += ftu;
		}
		else if (ftv > 0) {
			// Bottom face 
			float d = dev_sqrt(nx * nx + ny * ny);
			v = PI / 2 - dev_atan2f(d, (*dev_ak));
			u = PI / 2 - dev_atan2f(ny, -nx);
		}
		else {
			// Top face
			float d = dev_sqrt(nx * nx + ny * ny);
			v = -PI / 2.0 + dev_atan2f(d, (*dev_ak));
			u = -PI / 2.0 + dev_atan2f(-ny, nx);
		}
		u = u / (PI);
		v = v / (PI / 2.0);

		while (v < -1) {
			v += 2;
			u += 1;
		}
		while (v > 1) {
			v -= 2;
			u += 1;
		}

		while (u < -1) {
			u += 2;
		}
		while (u > 1) {
			u -= 2;
		}

		u = u / 2.0 + 0.5;
		v = v / 2.0 + 0.5;

		u = u * (de_img.cols - 1);
		v = v * (de_img.rows - 1);

		mapx(X, Y) = (float)u;
		mapy(X, Y) = (float)v;

		int u1 = floor(u);
		int v1 = floor(v);
		int u2 = u1 + 1;
		int v2 = v1 + 1;

		if (u1 >= 0 && v1 >= 0 && u2 < diff_area.cols && v2 < diff_area.rows)
		{
			//diff_area(int(v), int(u)) = (float)((*i) + 1);
			diff_area(int(v1), int(u1)) = (float)((*i) + 1);
			diff_area(int(v2), int(u2)) = (float)((*i) + 1);
			diff_area(int(v1), int(u2)) = (float)((*i) + 1);
			diff_area(int(v2), int(u1)) = (float)((*i) + 1);

			//Index_x(int(v), int(u)) = Y;
			Index_x(int(v1), int(u1)) = Y;
			Index_x(int(v2), int(u2)) = Y;
			Index_x(int(v1), int(u2)) = Y;
			Index_x(int(v2), int(u1)) = Y;

			//Index_y(int(v), int(u)) = X;
			Index_y(int(v1), int(u1)) = X;
			Index_y(int(v2), int(u2)) = X;
			Index_y(int(v1), int(u2)) = X;
			Index_y(int(v2), int(u1)) = X;
		}
	}
}

void createCubeMapFace(const cv::cuda::GpuMat & inputimg, const cv::cuda::GpuMat & cube_images, const cv::cuda::GpuMat & diff_area, const cv::cuda::GpuMat & Index_x, const cv::cuda::GpuMat & Index_y, cv::Mat outimg[])
{

	cv::cuda::GpuMat mapx(cube_height, cube_width, CV_32FC1);
	cv::cuda::GpuMat mapy(cube_height, cube_width, CV_32FC1);

	float an = sin(PI / 4.0);
	float ak = cos(PI / 4.0);

	cudaError_t Error;
	float* dev_an; float* dev_ak;
	Error = cudaMalloc((void**)&dev_an, sizeof(float));
	Error = cudaMemcpy(dev_an, &an, sizeof(float), cudaMemcpyHostToDevice);
	Error = cudaMalloc((void**)&dev_ak, sizeof(float));
	Error = cudaMemcpy(dev_ak, &ak, sizeof(float), cudaMemcpyHostToDevice);

	HOD_Point2f *gpu_faceTransform;

	HOD_Point2f faceTransform[6] = { HOD_Point2f(0.0, 0.0),		  //front
		HOD_Point2f(PI / 2.0, 0.0),  //right
		HOD_Point2f(PI, 0.0),		  //back
		HOD_Point2f(-PI / 2.0, 0.0), //left
		HOD_Point2f(0.0, -PI / 2.0), //top
		HOD_Point2f(0.0, PI / 2.0) };//bottom

	dim3 thread(32, 32);
	dim3 block((mapx.cols + thread.x - 1) / thread.x, (mapx.rows + thread.y - 1) / thread.y);

	Error = cudaMalloc((void**)&gpu_faceTransform, sizeof(HOD_Point2f) * 6);
	Error = cudaMemcpy(gpu_faceTransform, faceTransform, sizeof(HOD_Point2f) * 6, cudaMemcpyHostToDevice);
	int* dev_i;

	for (int i = 0; i < 6; i++)
	{
		Error = cudaMalloc((void**)&dev_i, sizeof(int));
		Error = cudaMemcpy(dev_i, &i, sizeof(int), cudaMemcpyHostToDevice);

		kernel7 << > > (inputimg, gpu_faceTransform, dev_i, dev_an, dev_ak, mapx, mapy, diff_area, Index_x, Index_y);
		cudaThreadSynchronize();
		//Mat cpu_mapx, cpu_mapy;
		//mapx.download(cpu_mapx);
		//mapy.download(cpu_mapy);
		//writeMapToFile(cpu_mapx, cpu_mapy, i);
		//Mat cpu_area;
		//diff_area.download(cpu_area);
		//diff_area_color(cpu_area);
		//Mat cpu_index_x, cpu_index_y;
		//Index_x.download(cpu_index_x);
		//Index_y.download(cpu_index_y);
		//writeEqu2cubeToFile(cpu_index_x, cpu_index_y);
		//Mat cpu_area;
		//diff_area.download(cpu_area);
		//writefaceindexToFile(cpu_area);

		cv::cuda::remap(inputimg, cube_images, mapx, mapy, CV_INTER_LINEAR, cv::BORDER_CONSTANT, cv::Scalar(0, 0, 0));

		cube_images.download(outimg[i]);
	}

	cudaFree(dev_an);
	cudaFree(dev_ak);
	cudaFree(dev_i);
	cudaFree(gpu_faceTransform);
}

打印图像outimg
writeCubeFace(cubeout);

void writeCubeFace(cv::Mat outimg[])
{
	for (int i = 0; i < 6; i++) {
		if (i == 0)
		{
			imwrite("./Image_out/segmentation/Cubemaps/front.jpg", outimg[i]);
		}
		else if (i == 1)
		{
			imwrite("./Image_out/segmentation/Cubemaps/right.jpg", outimg[i]);
		}
		else if (i == 2)
		{
			imwrite("./Image_out/segmentation/Cubemaps/back.jpg", outimg[i]);
		}
		else if (i == 3)
		{
			imwrite("./Image_out/segmentation/Cubemaps/left.jpg", outimg[i]);
		}
		else if (i == 4)
		{
			imwrite("./Image_out/segmentation/Cubemaps/top.jpg", outimg[i]);
		}
		else if (i == 5)
		{
			imwrite("./Image_out/segmentation/Cubemaps/bottom.jpg", outimg[i]);
		}
	}
}