计算机视觉-图像缩放算法-cuda实现

一.CUDA

CUDA是显卡厂商NVIDIA（英伟达）推出的运算平台，能够将数据数据复制到GPU，在GPU中进行计算，然后再返回给CPU端。CUDA将GPU称为设备侧或者Device，将CPU称为Host侧，按照下图进行划分，真正进行计算的是其中的Thread。

CUDA的极简使用流程
1.写核函数，是指在Device上运行的函数，用__global__ 标识。
2.在Device上申请空间，使用CUDA自带的cudaMalloc函数。
3.把要进行计算的数据拷贝到Device，也就是上一步申请到的空间。使用cudaMemcpy函数。
4.在Host中调用核函数，使用Device进行计算。
5.将计算的结果从Device拷贝回Host。

二.图像缩放算法

图像缩放这个事，所有操作图片的人都做过，但是程序是如何计算的？
首先，确定缩放目标图的大小，创建一张空白的目标图。
然后，从原图中为目标图的每一个像素点找到对应的值。
成熟的缩放算法有很多，在这里只介绍其中最简单的两个。
最近邻点插值法和双线性插值法

三.最近邻点插值法

1.算法

1.针对目标图中每一个像素点，按照缩放比例去原图中找对应的点，比值不为整数的就选目标图中最近的点。
2.给目标图取对应的值。

源图像与目标图像宽的比例为：

srcWidth / dstWidth

源图像与目标图像高的比例为：

srcHeight / dstHeight

由目标图坐标找原图像坐标就是：

srcX = dstX * (srcWidth / dstWidth)
srcY = dstY * (srcHeight / dstHeight)

2.代码

代码使用的是C++、CUDA

//文件名是test1.cu
#include<opencv2/opencv.hpp>
#include<iostream>
#include <string>
#include <stdio.h>
#include<vector>
using namespace cv;
using namespace std;

__global__ void resize1GPU(const unsigned char* src, int srcWidth, int srcHeight, unsigned char* dst, int dstWidth, int dstHeight)
{
	//核函数会在每个thread上运行，这里求的x、y是当前thread的坐标，同时也代表当前要处理的像素的坐标
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	int x = blockIdx.x * blockDim.x + threadIdx.x;	
	if (x >= dstWidth || y >= dstHeight) return;
	//以指针的形式操作图像，outPosition是指目标图像素在内存中的位置 
	int outPosition = y * dstWidth + x;		
	//求取对应原图的像素点，srcPosition是指原图像素在内存中的位置 
	int srcX = x * srcWidth / dstWidth;		//如果出现浮点数，这里就会向下取整，以此来表示最近邻
	int srcY = y * srcHeight / dstHeight;	//（如果不喜欢向下取整，也可以选择四舍五入）
	int srcPosition = srcY * srcWidth + srcX;		
	//为目标图像素赋值。RGB三通道，在内存中的位置是挨着的。
	dst[outPosition * 3 + 0] = src[srcPosition * 3 + 0];		
	dst[outPosition * 3 + 1] = src[srcPosition * 3 + 1];
	dst[outPosition * 3 + 2] = src[srcPosition * 3 + 2];
}

int main() {
	Mat src = imread("1.jpg");		//使用opencv
	int srcWidth = src.cols;
	int srcHeight = src.rows;
	int dstWidth = 256;					//目标图的大小
	int dstHeight = 256;

	unsigned char* devSrc;
	unsigned char* devDst;
	//在GPU上为两张图申请存储空间
	cudaMalloc((void**)&devSrc, srcWidth * srcHeight * 3 * sizeof(unsigned char));				
	cudaMalloc((void**)&devDst, dstWidth * dstHeight * 3 * sizeof(unsigned char));
	//把原图复制到GPU上，注意图片数据格式的变化
	cudaMemcpy(devSrc, (unsigned char*)(src.data), srcWidth * srcHeight * 3 * sizeof(unsigned char), cudaMemcpyHostToDevice);

	dim3 blocks((dstWidth + 15) / 16, (dstHeight + 15) / 16);
	dim3 threads(16, 16);
	//调用核函数，重点关注blocks与threads的设置，这样设置是为了让thread的坐标代表目标图像素的坐标
	resize1GPU <<<blocks, threads>>> (devSrc, srcWidth, srcHeight, devDst, dstWidth, dstHeight);
	//将处理完的目标图拷贝回来
	Mat dst(Size(dstWidth, dstHeight), CV_8UC3);
	cudaMemcpy(dst.data, devDst, dstWidth * dstHeight * 3 * sizeof(unsigned char), cudaMemcpyDeviceToHost);
	//使用opencv保存新图片
    vector<int> comprocession_params;
    comprocession_params.push_back(IMWRITE_PNG_COMPRESSION);
    comprocession_params.push_back(9);
    imwrite("nearest_neighbor.png", dst,comprocession_params);
	cudaFree(devSrc);
	cudaFree(devDst);
	return 0;
}

四.双线性插值法

1.算法

双线性插值的像素对应公式也是根据缩放比例
不一样的是，不再是找原图的1个点，而是找到原图的4个点，然后再求出目标图的像素值。
先思考这样一个问题：如下图，已知下图中所有点的坐标，并且已知Q12，Q22，Q11，Q21四点的值，假设所有点的值符合线性变化，那么如何求P点的值。

由线性假设有以下等式：

(R1-Q11) / (Q21-Q11) = (x-x1) / (x2-x1)
(R2-Q12) / (Q22-Q12) = (x-x1) / (x2-x1)
(R2-P) / (R2-R1) = (y2-y) / (y2-y1)

先由前两个等式求出R1，R2，做了一次线性插值
再由R1，R2求出P，又做了一次线性插值，这也是算法名字的由来。

最后求得
把这个求值算法转移到图像中，将坐标变成离散的情况。Q的四个点的坐标之间相差1，故进一步化简为：

P=(Q11)(x2- x)(y2-y) + (Q21)(x- x1)(y2- y) + (Q12)(x2- x)(- y1) +
(Q22)(x - x1)(y- y1)

1.由目标图中的像素，算出原图中的坐标（浮点数），一定能找到距离这个坐标最近的四个像素点。
2.用双线性插值计算出这点的值。

2.代码

//文件名是test2.cu
#include<opencv2/opencv.hpp>
#include<iostream>
#include <string>
#include <stdio.h>
#include<vector>
using namespace cv;
using namespace std;

__global__ void resize2GPU(const unsigned char* src, int srcWidth, int srcHeight, unsigned char* dst, int dstWidth, int dstHeight)
{
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	if (x >= dstWidth || y >= dstHeight) return;
	int dstOffset = y * dstWidth + x;	//目标图像素在内存中的位置 
	//根据缩放比例，计算在原图的坐标（浮点值）
	int srcXf = x * ((float)srcWidth / dstWidth);
	int srcYf = y * ((float)srcHeight / dstHeight);
	//向下取整，得到四个像素中，左上的像素坐标
	int srcX = (int)srcXf;
	int srcY = (int)srcYf;		
	//u就是上面算法中的x-x1，1-u就是x2-x
	int u = srcXf - srcX;		
	int v = srcYf - srcY;

	//P=(Q11)(x2- x)(y2-y) + (Q21)(x- x1)(y2- y) + (Q12)(x2- x)(- y1) + (Q22)(x - x1)(y- y1)
	dst[dstOffset] = 0;
	dst[dstOffset] += (1 - u) * (1 - v) * src[(srcY * srcWidth + srcX)];	
	dst[dstOffset] += (1 - u) * v * src[((srcY + 1) * srcWidth + srcX)];
	dst[dstOffset] += u * (1 - v) * src[(srcY * srcWidth + srcX + 1)];
	dst[dstOffset] += u * v * src[((srcY + 1) * srcWidth + srcX + 1)];	
	//(srcY+1)*srcWidth+srcX+1)是右下角的像素点在内存中的位置
}
//主函数和上一个算法代码一样，唯一区别就是，为了代码简单，把图片变成了灰度图
int main() {
	Mat src = imread("1.jpg", 0);
	int srcWidth = src.cols;
	int srcHeight = src.rows;
	int dstWidth = 512;
	int dstHeight = 512;

	unsigned char* devSrc;
	unsigned char* devDst;

	cudaMalloc((void**)&devSrc, srcWidth * srcHeight * sizeof(unsigned char));
	cudaMalloc((void**)&devDst, dstWidth * dstHeight * sizeof(unsigned char));
	cudaMemcpy(devSrc, src.data, srcWidth * srcHeight * sizeof(unsigned char), cudaMemcpyHostToDevice);

	dim3 blocks((dstWidth + 15) / 16, (dstHeight + 15) / 16);
	dim3 threads(16, 16);
	resize2GPU << <blocks, threads >> > (devSrc, srcWidth, srcHeight, devDst, dstWidth, dstHeight);

	Mat dst(Size(dstWidth, dstHeight), CV_8UC1);
	cudaMemcpy(dst.data, devDst, dstWidth * dstHeight * sizeof(unsigned char), cudaMemcpyDeviceToHost);

	vector<int> comprocession_params;
    comprocession_params.push_back(IMWRITE_PNG_COMPRESSION);
    comprocession_params.push_back(9);
    imwrite("Bilinear_Interpolation.png", dst,comprocession_params);
	cudaFree(devSrc);
	cudaFree(devDst);
	return 0;
}

3.遗留问题

1.当出现坐标在边界时，会出现(i+1,j)、(i,j+1)、(i+1,j+1)的点在实际的图中不存在？此时可以用(i,j)、(i-1,j)、(i,j-1)、(i-1,j-1)。
2.当按缩放比计算坐标后，与原图像中的点发生重合，怎么选取剩下3个点？无论我们怎么选取，其实其余3点的权重都是0，所以都不影响。
3.目标图像的原点(0, 0)点和原始图像的原点肯定是重合的，但是目标图的中心点不一定与原图的中心点重合。整体图像会发生偏移。把浮点坐标的计算方法改为：

4.我的电脑没有英伟达显卡，租了一个云计算平台来运行.cu格式的代码。
使用以下两个命令行
编译：nvcc test1.cu pkg-config opencv --libs --cflags opencv -o test1
运行：./test1

CUDA-图像缩放