http://peghoty.blog.163.com/blog/static/493464092013016113254852/
http://blog.csdn.net/augusdi/article/details/12833235
CUDA存储器模型：http://blog.csdn.net/endlch/article/details/44538801
CUDA限定符：http://blog.csdn.net/shouhuxianjian/article/details/42427285
思想即是将内存数据拷贝到显存，在显存上执行并行运算，将结果数据从显存拷贝回内存。
CUDA内有thrust库，类似于C++ stl库。
===========以下是原文=========
挖坑待填。
以上是本机CUDA参数。
需要了解的概念：线程束(wrap)，共享内存，常量内存，纹理内存(?，图形学相关，略)，流，原子操作。
寄存器
寄存器是GPU片上高速缓存， 执行单元可以以极低的延迟访问寄存器。寄存器的基本单元式寄存器文件，每个寄存器文件大小为32bit。局部存储器对于每个线程，局部存储器也是私有的。如果寄存器被消耗完。数据将被存储在局部存储器中。如果每个线程使用了过多的寄存器，或声明了大型结构体或数据，或者编译器无法确定数据的大小，线程的私有数据就有可能被分配到local memory中，一个线程的输入和中间变量将被保存在寄存器或者是局部存储器中。局部存储器中的数据被保存在显存中，而不是片上的寄存器或者缓存中，因此对local memory的访问速度很慢。
共享存储器
共享存储器(share memeory)也是GPU片内缓存存储器。它是一块可以被同一block中的所有线程访问的可读存储器。
使用关键字share添加到变量的声明中，这将使这个变量驻留在共享内存中。cuda c编译器对共享内存中的变量与普通变量将采取不同的处理方式。
对于在GPU上启动的每个线程块，cuda c编译器都将创建该变量的一个副本，线程块中的每一个线程都共享这块内存，但这个线程却无法看到也不能修改其他线程块的变量的副本。这就实现了一种非常好的方式，使得一个线程块中的多个线程能够在计算上进行通信和协作，而且，共享内存缓冲区驻留在物理GPU上，而不是驻留在GPU之外的系统内存中。
常量内存
__constant__将把变量的访问限制为只读。
在接受了这种限制之后，我们希望得到某种回报，与全局内存中读数据相比，从常量内存中读取相同的数据可以节约内存的带宽，主要有两个原因：
-对常量内存的单次读操作可以广播到其他的“领进”线程，这将节约15次读取操作。
-常量内存的数据缓存起来，因此对相同地址的连续读取操作将不会产生额外的内存通信量。
“邻近”是指半个warp中的线程。当处理常量内存时。nvidia硬件将把单次内存读取操作广播到每个半线程束。在半线程束中包含了16个线程，即线程束中数量的一半。
如果在半线程束中的每一个线程访问相同的常量内存地址。那么GPU只会发生一次读操作事件并在随后将数据广播到每个线程。
如果从常量内存中读取大量的数据，那么这种方式产生的内存流量只是全局内存时的1/16.然而，当使用常量内存时也可能产生负面影响。
如果半线程束的所有16个线程需要访问常量内存中不同的数据，那么这个16次读取操作会被串行化，从而需要16倍的时间来发出请求。
但如果从全局内存中读取，那么这些请求会同时发出。在这种情况下，从常量内存读取就慢于从全局内存中读取。
全局存储器
全局存储器(global memeory)位于显存（占据了大部分的显存)。
整个网格中的任意线程都能读写全局存储器的任意位置。在目前的架构中，全局存储器没有缓存。
=======补======
SM与线程束
通常，线程块的数量为GPU*处理器数量的2倍时，将达到最优性能。
GPU拥有数百个核，其中，SM代表多流处理器，即计算核心，而每个SM又包含8个标准流处理器SP，以及其他。
隶属于同一个SM的8个SP共用同一套取指与发射单元，共用一块共享存储器。
kernel以block为单位执行。
一个block必须被分配到同一块SM中，block的每个thread会发送到一个SP上执行，但一个SM中同一时刻可以有多个活动线程块在等待执行。
同一个block的thread开始于相同的指令地址，理论上能够按不同的分支执行，但实际上由于8个SP共用一套取值与发射单元，同一warp的线程执行的指令是相同的。
如果一个warp的线程跳转如分支语句的同一分支，那么实际执行时间就是这个分支执行时间；
否则，SM需要把每一个分支的指令发射到每个SP，执行时间是执行多个分支的所用时间之和。
故CUDA程序尽量避免分支，尽量warp内不分支。
线程束(warp)：一个线程束由32个连续的线程组成。（简单地说，warp包含32个线程是因为每发射一条warp指令，SM中的8个SP就会将这条指令执行4遍）。warp才是真正的执行单位。
杂
原子操作。同时对global内存写操作，可分批进行，改成先线程块对shared内存写操作，结束后shared内存写入global内存。
__syncthreads()实现了线程块内的线程同步，当任意线程运行到BAR标记处后，暂停运行，直到整个block中所有的thread都运行到BAR标记处后才继续执行。
__syncthreads()勿置于分支语句内。
流：名义上多个流，实际上可能就是kenel（GPU运算）和copy（CPU与GPU间数据传送）两个流。每个流都是一个队列，事件被push入队列中等待执行。
for循环的任务切分的时候，有两种方式划分任务。
1.划分成k段，每段由某个线程执行。
2.按模k同余进行划分，for循环每次的递增量为块大小。
一般第2种方式更优，因为是并行执行，故第二种方式保证每次执行的的时候，不同线程访问的数据位置相邻。
并行算法
归约运算: 每次折半。以求和为例，第一次前1/2 + 后1/2；第二次 前1/4 + 后1/4 .。。
int i = blockDim.x/;
while(i != ) {
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
i /= ;
}
if (cacheIndex == )
c[blockIdx.x] = cache[];
更好的优化算法：循环展开等
前缀和运算(Scan):
for(d = ; (1 << d) < n; d++)
for all k in parallel
if( k >= (1 << d) )
x[out][k] = x[in][k – (1 << (d-1))] + x[in][k]
else
x[out][k] = x[in][k]
或
for(d = ; ( << d) < n; d++)
for all k in parallel
tmp = x[k]
if( k >= ( << d) )
tmp = x[k – ( << (d-))] + x[k]
__syncthreads();//同步
x[k] = tmp
以上两算法运行所需空间至少是原空间的两倍，思想为倍增思想。
还有更高效的Scan算法。
for d:= to log2(n-) do
for k from to n- by ^(d+) in parallel do
x[k+^(d+)-]:=x[k+^(d+)-] + x[k+^(d+)-] x[n-]:=
for d:=log2(n-) downto do
for k from to n- by ^(d+) in parallel do
t:=x[k+^d-]
x[k+^d-]:=x[k+^(d+)-]
x[k+^(d+)-]:=t+x[k+^(d+)-]
书上还有更高效的scan_best_kernel算法，略。
排序算法:
基于比较的排序：排序网络
基于非比较的排序：并行基数排序。前置技能：Scan。
并行基数排序算法：
按二进制位考虑。
以00101101为例。排完序后应当是12473568。
二进制翻转:
统计前缀和:
如果当前位是0，则写入对应位置。
第1个数写入首位置，第2个数写入第二个位置，第4个数写入第三个位置，第7个数写入第四个位置。
再对当前位是1的进行写入，位置下标 + (0的个数)。
矩阵乘法优化:
矩阵运算A*B = C， A为m*p的矩阵，B为p*n的矩阵。
优化1:
将C分块，每个线程块处理C中的某一块，例如d*d的一小块。
那么每个线程块需要完成d*p的矩阵与p*d的矩阵相乘的运算。
为了高效访存，每个线程块再对d*p和p*d的矩阵的p进行划分，看成多个矩阵块相乘后累加。
每个小块为d*q和q*d的大小，开在shared memory内，节约了大量global memory带宽。
（虽然循环次数会增加，但访存效率得到了高效提升）
优化2：
利用寄存器资源优化，效率更高，但略为繁琐。
矩阵转置优化:
无优化：拷贝至GPU内存，置换后拷贝回CPU内存。缺点：输入时每个block按行读入，满足合并访问条件；输出时数据间隔过大，不满足合并访问条件。
优化1：
分块，每个块是一个小方阵矩阵，如16*16。
输入时，每个线程块操作一个16*16方阵，通过shared memory完成16*16小方阵转置。
之后将大矩阵按块转置输出至global memory，每个线程块内无需再转置，满足合并访问条件。
shared memory数组大小设置成16*17而不是16*16，这样每行中处于同一列的数据就会被存储在不同的shared memory bank中，避免了bank conflict。
优化2：
上述无优化与优化1均存在分区冲突问题。优化2算法进行了for循环操作，暂未深入研究。
CUDA程序优化
grid和block的维度设计
grid的尺寸大一点较好。
为了有效利用执行单元，每个block的线程数应当是32的整数倍，最好让线程数量保持在64 ~ 256之间。
block维度和每个维度上的尺寸的主要作用是避免做整数除法和求模运算。实际中视具体情况而定。
如果问题规模对划分方式不敏感，应该让blockDim.x为16或16的整数倍，提高访问global memory和shared memory的效率。
存储器访问优化
灵活运用各种存储器的特性，实现最大可用带宽。
指令流优化
等
CUDA作业
作业1 简单CUDA应用，矩阵乘法
#include <bits/stdc++.h>
//#include "cuda_runtime.h"
//#include "device_launch_parameters.h" using namespace std;
#define N 2000 const int block = <<;
const int thread = <<; long long a[N][N];
long long b[N][N];
long long c[N][N];
void init() {
for(int i = ; i < N; i++)
for(int j = ; j < N; j++)
a[i][j] = i*i+j, b[i][j] = i+j*j, c[i][j] = ;
} __global__ void init_cuda(long long *c) {
int id = blockIdx.x*blockDim.x+threadIdx.x;
if(id < N*N) c[id] = ;
} __global__ void mul_cuda(long long *a, long long *b, long long *c) {
int id = blockIdx.x*blockDim.x+threadIdx.x;
if(id < N*N) {
int row = id/N, col = id-row*N;
for(int k = ; k < N; k++)
c[id] += a[row*N+k]*b[k*N+col];
}
} int main(int argc, char** argv) {
int cstart = clock();
init();
if(argv[][] == '') {
puts("not cuda");
for(int i = ; i < N; i++)
for(int j = ; j < N; j++)
for(int k = ; k < N; k++)
c[i][k] += a[i][j]*b[j][k];
}
else {
puts("cuda");
long long *dev_a, *dev_b, *dev_c;
cudaMalloc( (void**)&dev_a, sizeof a );
cudaMemcpy(dev_a, a, sizeof a, cudaMemcpyHostToDevice); cudaMalloc( (void**)&dev_b, sizeof b );
cudaMemcpy(dev_b, b, sizeof b, cudaMemcpyHostToDevice); cudaMalloc( (void**)&dev_c, sizeof c ); init_cuda<<<block, thread>>>(dev_c);
mul_cuda<<<block, thread>>>(dev_a, dev_b, dev_c); cudaMemcpy(c, dev_c, sizeof c, cudaMemcpyDeviceToHost);
cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_c);
}
printf("%lld, ", c[][]);
printf("time: %d\n", int(clock()-cstart));
return ;
}
作业2 卷积操作，常量内存
//compile command: nvcc cv.cu `pkg-config --cflags --libs opencv` -std=c++11
//execute command1: ./a.out CC.jpg 3
//execute command2: ./a.out CC.jpg 5
#include <bits/stdc++.h> #include <opencv2/opencv.hpp>
//#include <opencv2/gpu/gpu.hpp>
using namespace cv; Mat G3 = (Mat_<int>(, ) << -, , -,
, , ,
-, , -); Mat G5 = (Mat_<int>(, ) << , -, , -, ,
-, , -, ,-,
, -, ,-, ,
-, , -, ,-,
, -, ,-, ); void CPU_Sharpen(const Mat& myImage, Mat& Result, int ca){
CV_Assert(myImage.depth() == CV_8U); // accept only uchar images int begin = clock();
Result.create(myImage.size(), myImage.type());
const int nChannels = myImage.channels();
Mat &G = ca == ? G3: G5;
int half = G.rows >> ; for(int row = half; row < myImage.rows-half; ++row) {
uchar* output = Result.ptr<uchar>(row);
for(int col = half*nChannels; col < nChannels * (myImage.cols - half); ++col) {
int tmp = ;
for(int i = ; i < G.rows; ++i)
for(int j = ; j < G.cols; ++j)
tmp += G.at<int>(i, j)*( *(myImage.ptr<uchar>(row-half+i)+(col-half*nChannels+j*nChannels) ) );
*output++ = saturate_cast<uchar>(tmp);
}
}
for(int i = ; i < half; i++) {
Result.row(i).setTo(Scalar());
Result.row(Result.rows - - i).setTo(Scalar());
Result.col(i).setTo(Scalar());
Result.col(Result.cols - - i).setTo(Scalar());
}
printf("Time used %.3fms\n", ((int)clock()-begin)*1000.0/CLOCKS_PER_SEC);
} /********************************************/ __constant__ int con_G3[][] = {
{-, , -},
{ , , },
{-, , -}
};
__constant__ int con_G5[][] = {
{, -, ,-, },
{-, , -, ,-},
{, -, ,-, },
{-, , -, ,-},
{, -, ,-, }
}; __global__ void init_cuda(uchar *c, int col_num) {
int col_id = blockIdx.x, row_id = threadIdx.x;
int now = (row_id*col_num+col_id)*;
c[now] = c[now+] = c[now+] = ;
} //GPU,start from c, num * sizeof(uchar)
__global__ void test(uchar *c, int *sum, int num) {
int x = ;
for(int i = ; i < num; i++)
x += c[i];
*sum = x;
} __global__ void con_cuda(uchar *s, uchar *t, int ca, int row_num, int col_num) {
int col_id = blockIdx.x-, row_id = threadIdx.x-;
const int half = ca >> ;
if(row_id >= half && row_id < row_num-half && col_id >= half && col_id < col_num-half) {
const int* con_mat = ca == ? con_G3[]: con_G5[];
int res[] = {, , };
for(int i = ; i < ca; i++)
for(int j = ; j < ca; j++) {
//s[row_num][col_num][3];
int pos = (row_id-half+i)*col_num*+(col_id-half+j)*;
res[] += con_mat[i*ca+j]*s[pos];
res[] += con_mat[i*ca+j]*s[pos+];
res[] += con_mat[i*ca+j]*s[pos+];
}
res[] = res[] < ? : (res[] > ? : res[]);
res[] = res[] < ? : (res[] > ? : res[]);
res[] = res[] < ? : (res[] > ? : res[]);
int pos = row_id*col_num*+col_id*;
t[pos] = res[],
t[pos+] = res[],
t[pos+] = res[];
}
} /*******************************************/ void HANDLE(cudaError x) {
if(x != cudaSuccess) {
puts("error!");
exit();
}
} int main(int argc, char** argv ) {
if ( argc < ) {
printf("usage: a.out <Image_Path> <size of Mat>\n");
return -;
} Mat src_img = imread(argv[], ), ans_CPU, ans_GPU;
int ca = argv[][]-'';
printf("%d %d\n", src_img.rows, src_img.cols); /**********************************************************************************/ printf("Run on CPU!\n");
CPU_Sharpen(src_img, ans_CPU, ca);
std::string s = std::string("CC")+std::to_string(ca)+std::string("_With_CPU.jpg");
imwrite(s, ans_CPU);
imshow("after operation", ans_CPU);
waitKey(); /**********************************************************************************/ printf("Run on GPU!\n");
int begin = clock();
uchar *dev_src, *dev_result;
int seg = src_img.cols*src_img.channels();
HANDLE(cudaMalloc( (void**)&dev_src, src_img.rows*seg*sizeof(uchar)));
HANDLE(cudaMalloc( (void**)&dev_result, src_img.rows*seg*sizeof(uchar)));
/*Memcpy to dev_src*/
for(int i = ; i < src_img.rows; ++i)
HANDLE(cudaMemcpy(dev_src+i*seg*sizeof(uchar), src_img.ptr<uchar>(i), sizeof(uchar)*seg, cudaMemcpyHostToDevice));
/*Init for dev_result*/
init_cuda<<<src_img.cols, src_img.rows>>>(dev_result, src_img.cols);
/*Do convolution*/
con_cuda<<<src_img.cols, src_img.rows>>>(dev_src, dev_result, ca, src_img.rows, src_img.cols); ans_GPU.create(src_img.size(), src_img.type());
/*Memcpy to host*/
for(int i = ; i < ans_GPU.rows; ++i)
cudaMemcpy(ans_GPU.ptr<uchar>(i), dev_result+i*seg*sizeof(uchar), sizeof(uchar)*seg, cudaMemcpyDeviceToHost); for(int i = ; i < (ca >> ); i++) {
ans_GPU.row(i).setTo(Scalar());
ans_GPU.row(ans_GPU.rows - - i).setTo(Scalar());
ans_GPU.col(i).setTo(Scalar());
ans_GPU.col(ans_GPU.cols - - i).setTo(Scalar());
}
/*Free*/
cudaFree(dev_src);
cudaFree(dev_result);
printf("Time used %.3fms\n", ((int)clock()-begin)*1000.0/CLOCKS_PER_SEC);
imshow("after operation", ans_GPU);
waitKey();
return ;
}
作业3 卷积操作，流
//compile command: nvcc cv.cu `pkg-config --cflags --libs opencv` -std=c++11
//execute command1: ./a.out 1.jpg 3
//execute command2: ./a.out 1.jpg 5
#include <bits/stdc++.h> #include <opencv2/opencv.hpp>
//#include <opencv2/gpu/gpu.hpp>
using namespace cv; /********************************************/ __constant__ int con_G3[][] = {
{-, , -},
{ , , },
{-, , -}
};
__constant__ int con_G5[][] = {
{, -, ,-, },
{-, , -, ,-},
{, -, ,-, },
{-, , -, ,-},
{, -, ,-, }
}; __global__ void init_cuda(uchar *c, int col_num) {
int col_id = blockIdx.x, row_id = threadIdx.x;
int now = (row_id*col_num+col_id)*;
c[now] = c[now+] = c[now+] = ;
} //GPU,start from c, num * sizeof(uchar)
__global__ void test(uchar *c, int *sum, int num) {
int x = ;
for(int i = ; i < num; i++)
x += c[i];
*sum = x;
} __global__ void con_cuda(uchar *s, uchar *t, int ca, int row_num, int col_num) {
int col_id = blockIdx.x-, row_id = threadIdx.x-;
const int half = ca >> ;
if(row_id >= half && row_id < row_num-half && col_id >= half && col_id < col_num-half) {
const int* con_mat = ca == ? con_G3[]: con_G5[];
int res[] = {, , };
for(int i = ; i < ca; i++)
for(int j = ; j < ca; j++) {
//s[row_num][col_num][3];
int pos = (row_id-half+i)*col_num*+(col_id-half+j)*;
res[] += con_mat[i*ca+j]*s[pos];
res[] += con_mat[i*ca+j]*s[pos+];
res[] += con_mat[i*ca+j]*s[pos+];
}
res[] = res[] < ? : (res[] > ? : res[]);
res[] = res[] < ? : (res[] > ? : res[]);
res[] = res[] < ? : (res[] > ? : res[]);
int pos = row_id*col_num*+col_id*;
t[pos] = res[],
t[pos+] = res[],
t[pos+] = res[];
}
} /*******************************************/ void HANDLE_ERROR(cudaError x) {
if(x != cudaSuccess) {
puts("error!");
exit();
}
} int main(int argc, char** argv ) {
if ( argc < ) {
printf("usage: a.out <Image_Path> <size of Mat>\n");
return -;
} Mat src_img = imread(argv[], ), ans_CPU, ans_GPU;
int ca = argv[][]-'';
printf("%d %d\n", src_img.rows, src_img.cols);
/**********************************************************************************/ printf("Run on GPU!\n");
int begin = clock();
/**********************************************************************************/ uchar *dev_src, *dev_result;
int seg = src_img.cols*src_img.channels();
HANDLE_ERROR(cudaMalloc( (void**)&dev_src, src_img.rows*seg*sizeof(uchar)));
HANDLE_ERROR(cudaMalloc( (void**)&dev_result, src_img.rows*seg*sizeof(uchar))); cudaStream_t stream0, stream1;
HANDLE_ERROR( cudaStreamCreate( &stream0 ) );
HANDLE_ERROR( cudaStreamCreate( &stream1 ) );
for(int i = ; i < src_img.rows; ++i)
HANDLE_ERROR( cudaMemcpyAsync(dev_src+i*seg*sizeof(uchar), src_img.ptr<uchar>(i), sizeof(uchar)*seg, cudaMemcpyHostToDevice, stream0) ); init_cuda<<<src_img.cols, src_img.rows, , stream0>>>(dev_result, src_img.cols);
con_cuda<<<src_img.cols, src_img.rows, , stream0>>>(dev_src, dev_result, ca, src_img.rows, src_img.cols); ans_GPU.create(src_img.size(), src_img.type()); for(int i = ; i < ans_GPU.rows; ++i)
HANDLE_ERROR( cudaMemcpyAsync(ans_GPU.ptr<uchar>(i), dev_result+i*seg*sizeof(uchar), sizeof(uchar)*seg, cudaMemcpyDeviceToHost, stream0) ); for(int i = ; i < (ca >> ); i++) {
ans_GPU.row(i).setTo(Scalar());
ans_GPU.row(ans_GPU.rows - - i).setTo(Scalar());
ans_GPU.col(i).setTo(Scalar());
ans_GPU.col(ans_GPU.cols - - i).setTo(Scalar());
} HANDLE_ERROR( cudaStreamSynchronize( stream0 ) );
HANDLE_ERROR( cudaStreamSynchronize( stream1 ) ); HANDLE_ERROR( cudaStreamDestroy( stream0 ) );
HANDLE_ERROR( cudaStreamDestroy( stream1 ) ); cudaFree(dev_src);
cudaFree(dev_result);
/********************without stream*********************/ // uchar *dev_src, *dev_result;
// int seg = src_img.cols*src_img.channels();
// HANDLE_ERROR(cudaMalloc( (void**)&dev_src, src_img.rows*seg*sizeof(uchar)));
// HANDLE_ERROR(cudaMalloc( (void**)&dev_result, src_img.rows*seg*sizeof(uchar)));
// /*Memcpy to dev_src*/
// for(int i = 0; i < src_img.rows; ++i)
// HANDLE_ERROR(cudaMemcpy(dev_src+i*seg*sizeof(uchar), src_img.ptr<uchar>(i), sizeof(uchar)*seg, cudaMemcpyHostToDevice));
// /*Init for dev_result*/
// init_cuda<<<src_img.cols, src_img.rows>>>(dev_result, src_img.cols);
// /*Do convolution*/
// con_cuda<<<src_img.cols, src_img.rows>>>(dev_src, dev_result, ca, src_img.rows, src_img.cols);
//
// ans_GPU.create(src_img.size(), src_img.type());
// /*Memcpy to host*/
// for(int i = 0; i < ans_GPU.rows; ++i)
// HANDLE_ERROR( cudaMemcpy(ans_GPU.ptr<uchar>(i), dev_result+i*seg*sizeof(uchar), sizeof(uchar)*seg, cudaMemcpyDeviceToHost) );
//
// for(int i = 0; i < (ca >> 1); i++) {
// ans_GPU.row(i).setTo(Scalar(140));
// ans_GPU.row(ans_GPU.rows - 1 - i).setTo(Scalar(140));
// ans_GPU.col(i).setTo(Scalar(140));
// ans_GPU.col(ans_GPU.cols - 1 - i).setTo(Scalar(140));
// } /*Free*/
// cudaFree(dev_src);
// cudaFree(dev_result);
printf("Time used %.3fms\n", ((int)clock()-begin)*1000.0/CLOCKS_PER_SEC);
imshow("after operation", ans_GPU);
imwrite("Tigerwith5.jpg", ans_GPU);
waitKey();
return ;
}
final project 图墙 + 图片融合
//compile command: nvcc final.cu `pkg-config --cflags --libs opencv` -std=c++11
//execute command1: ./a.out CC.jpg 3
//execute command2: ./a.out CC.jpg 5
#include <bits/stdc++.h> #include <opencv2/opencv.hpp>
//#include <opencv2/gpu/gpu.hpp>
using namespace cv; __global__ void init_cuda(uchar *c, int col_num) {
int col_id = blockIdx.x, row_id = threadIdx.x;
int now = (row_id*col_num+col_id)*;
c[now] = c[now+] = c[now+] = ;
} //GPU,start from c, num * sizeof(uchar)
__global__ void test(uchar *c, int *sum, int num) {
int x = ;
for(int i = ; i < num; i++)
x += c[i];
*sum = x;
} __global__ void solve(uchar *s, uchar *t, uchar a, uchar b, uchar c, int R, int C) {
int x = blockIdx.x-, y = threadIdx.x-;
if(x < R&&y < C) {
int g = *(x*C+y);
if(t[g] == a&&t[g+] == b&&t[g+] == c) {
t[g] = s[g];
t[g+] = s[g+];
t[g+] = s[g+];
}
else {
t[g] = 0.35*s[g] +0.65*t[g];
t[g+] = 0.35*s[g+]+0.65*t[g+];
t[g+] = 0.35*s[g+]+0.65*t[g+];
}
}
} __global__ void zoom(uchar *s, uchar *t, int R, int C, int r, int c) {
int x = blockIdx.x-, y = threadIdx.x-;
if(x <= r && y <= c) {
int row = x/(float)r*R, col = y/(float)c*C;
//t[x][y] = s[row][col];
t[ (x*c+y)* ] = s[ (row*C+col)* ];
t[ (x*c+y)*+ ] = s[ (row*C+col)*+ ];
t[ (x*c+y)*+ ] = s[ (row*C+col)*+ ];
}
} /*******************************************/ void HANDLE(cudaError x) {
if(x != cudaSuccess) {
puts("error!");
exit();
}
} const int N = ; int main(int argc, char** argv ) {
Mat src_img[N], dst_img[N], ret1, ret2;
int width = , height = ;
for(int i = ; i < N; i++) {
src_img[i] = imread(std::to_string(i+)+std::string(".jpg"), );
int r = src_img[i].rows, c = src_img[i].cols;
if(height > r) height = r;
if(width > c) width = c;
}
height *= 0.4;
width *= 0.17; //resize
int begin = clock();
for(int i = ; i < N; i++) {
dst_img[i].create(Size(height, width), src_img[i].type());
resize(src_img[i], dst_img[i], Size(height, width));
}
printf("Time used in resizing is%.3fms\n", ((int)clock()-begin)*1000.0/CLOCKS_PER_SEC); int sq = sqrt(N+0.5);
//std::cout << width << ' ' << height << std::endl;
ret1.create(Size(height*sq, width*sq), src_img[].type());
//std::cout << ret1.rows << ' ' << ret1.cols << std::endl; //merge
for(int i = ; i < sq; i++)
for(int j = ; j < sq; j++){
for(int r = ; r < width; r++)
memcpy(ret1.ptr<uchar>(i*width+r)+j*height*, dst_img[i*sq+j].ptr<uchar>(r), height*);
} Mat ret = imread("0.jpg", );
resize(ret, ret2, Size(height*, width*));
//std::cout << ret2.rows << ' ' << ret2.cols << std::endl;
//imshow("", ret2);
//waitKey(); uchar a = *ret2.ptr<uchar>(), b = *(ret2.ptr<uchar>()+), c = *(ret2.ptr<uchar>()+);
if(ret1.rows != ret2.rows || ret1.cols != ret2.cols) puts("gg");
int R = ret2.rows, C = ret2.cols;
std::cout << R << ' ' << C << std::endl;
//CPU
begin = clock();
for(int i = ; i < R; i++) {
uchar *p1 = ret1.ptr<uchar>(i), *p2 = ret2.ptr<uchar>(i);
bool tag = true;
double x = ;
for(int j = ; j < C; j++) {
if(*(p2+j*) == a&&*(p2+j*+) == b&&*(p2+j*+) == c) {
x = ;
*(p2+j*) = *(p1+j*);
*(p2+j*+) = *(p1+j*+);
*(p2+j*+) = *(p1+j*+);
continue ;
} if(*(p2+j*+) == a&&*(p2+j*+) == b&&*(p2+j*+) == c) {
x = ;
*(p2+j*) = *(p1+j*);
*(p2+j*+) = *(p1+j*+);
*(p2+j*+) = *(p1+j*+);
continue ;
} x = tag? x+0.06: x-0.06;
if(x > ) tag = false;
if(x < 0.6) tag = true;
*(p2+j*) = (-x)*(*(p1+j*)) +x*(*(p2+j*));
*(p2+j*+) = (-x)*(*(p1+j*+))+x*(*(p2+j*+));
*(p2+j*+) = (-x)*(*(p1+j*+))+x*(*(p2+j*+));
}
}
printf("Time used in CPU is %.3fms\n", ((int)clock()-begin)*1000.0/CLOCKS_PER_SEC); imshow("", ret2);
waitKey();
imwrite("final.jpg", ret2);
//GPU
begin = clock();
uchar *dev_src, *dev_result;
HANDLE(cudaMalloc( (void**)&dev_src, R*C**sizeof(uchar)));
HANDLE(cudaMalloc( (void**)&dev_result, R*C**sizeof(uchar)));
for(int i = ; i < R; ++i)
cudaMemcpy(dev_src+i*C**sizeof(uchar), ret1.ptr<uchar>(i), C**sizeof(uchar), cudaMemcpyHostToDevice);
for(int i = ; i < R; ++i)
cudaMemcpy(dev_result+i*C**sizeof(uchar), ret2.ptr<uchar>(i), C**sizeof(uchar), cudaMemcpyHostToDevice);
solve<<<R, C>>>(dev_src, dev_result, a, b, c, R, C);
for(int i = ; i < R; ++i)
cudaMemcpy(ret2.ptr<uchar>(i), dev_result+i*C**sizeof(uchar), C**sizeof(uchar), cudaMemcpyDeviceToHost);
printf("Time used in GPU is %.3fms\n", ((int)clock()-begin)*1000.0/CLOCKS_PER_SEC); imshow("", ret2);
waitKey();
cudaFree(dev_src);
cudaFree(dev_result);
imwrite("final2.jpg", ret2); /**********************************************************************************/ printf("Run on CPU!\n");
CPU_Sharpen(src_img, ans_CPU, ca);
std::string s = std::string("IMG")+std::to_string(ca)+std::string("_With_CPU.jpg");
imwrite(s, ans_CPU);
imshow("after operation", ans_CPU);
waitKey(); /**********************************************************************************/ printf("Run on GPU!\n");
int begin = clock();
uchar *dev_src, *dev_result;
int seg = src_img.cols*src_img.channels();
HANDLE(cudaMalloc( (void**)&dev_src, src_img.rows*seg*sizeof(uchar)));
HANDLE(cudaMalloc( (void**)&dev_result, src_img.rows*seg*sizeof(uchar))); /*Memcpy to dev_src*/ for(int i = ; i < src_img.rows; ++i)
HANDLE(cudaMemcpy(dev_src+i*seg*sizeof(uchar), src_img.ptr<uchar>(i), sizeof(uchar)*seg, cudaMemcpyHostToDevice)); init_cuda<<<src_img.cols, src_img.rows>>>(dev_result, src_img.cols); con_cuda<<<src_img.cols, src_img.rows>>>(dev_src, dev_result, ca, src_img.rows, src_img.cols); ans_GPU.create(src_img.size(), src_img.type()); for(int i = ; i < ans_GPU.rows; ++i)
cudaMemcpy(ans_GPU.ptr<uchar>(i), dev_result+i*seg*sizeof(uchar), sizeof(uchar)*seg, cudaMemcpyDeviceToHost); for(int i = ; i < (ca >> ); i++) {
ans_GPU.row(i).setTo(Scalar());
ans_GPU.row(ans_GPU.rows - - i).setTo(Scalar());
ans_GPU.col(i).setTo(Scalar());
ans_GPU.col(ans_GPU.cols - - i).setTo(Scalar());
} cudaFree(dev_src);
cudaFree(dev_result);
printf("Time used %.3fms\n", ((int)clock()-begin)*1000.0/CLOCKS_PER_SEC);
imshow("after operation", ans_GPU);
waitKey(); return ;
}
emmmm，可能有bug，有冗余代码。
效果：
