1、并行算法介绍

Reduce:

对于一般可以加法运算，如进行8个数相加：

串行计算

由于对于符合交换律，其运算顺序可以改变，可以使用并行计算，如下:

并行计算

对比串行运算，与并行运算的速度：

运算量	8	...	N
并行Step	7	...	N-1
串行Step	3	...	log2(N)

对于并行计算的实现，可参照udacity的课程GitHub:reduce.cu

Scan

Scan 可称为prefix sum，即求前N项（或N-1项）的和，如下：

prefix sum
串行计算流程如上求和无太大差异，而并行计算起来相当讲技巧，以N=8举例:
黑\红实线：加法 黑虚线：赋值
详细可阅读：Parallel prefix sum (scan) with CUDA

2、图像的直方图均衡

假设存在8个灰度级的图像如下图，那么它实际可见的范围为1~4

1	2	3	4

1	2	3	4

1	2	3	4

1	2	3	4

频数图如下：

image.png histogram_equalize.png

即原来的灰度级为1映射后为0，2映射后为2，3映射后为3，4映射后为5
均衡后新灰度频数表如下：

灰度级	0	1	2	3	4	5	6	7
频数	4	0	4	4	0	4	0	0

image.png

映射结果：

0	2	3	5

0	2	3	5

0	2	3	5

0	2	3	5

总的来说没有整个过程没有改变像素点灰度值点之间大小、数量关系，或者说将可视区间拉伸，让色彩对比度更加明显。

3、并行图像直方图均衡化

1)统计每个灰度级像素的数目

__global__ void  histogram_init(char *src,unsigned int *out)
//针对512*512大小灰度图 开启512*512线程，计算256个灰度级相应频数

2)应用Parallel Scan(prefix sum) 计算出映射关系

__global__ void scan_downsweep(unsigned int *in,unsigned int *out,int n)
//输入256点频数表，输出均衡化后256对应得映射值

3)将原图按照映射关系映射到结果图

__global__ void equal_remap(char *img,char *out,unsigned int *idetity)
//根据映射表，重新映射图像

程序效果图如下：

histogram_result

源代码：

#include "cuda.h"
#include "cuda_runtime.h"
#include "cuda_runtime_api.h"
#include "device_functions.h"

#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/opencv.hpp"
#include "stdio.h"
#include <vector>
using namespace std;
using namespace cv;
#define IMAGE_DIR "/home/dzqiu/Documents/image/chaozhou.JPG"

__global__ void  histogram_init(char *src,unsigned int *out)
{
    int x = threadIdx.x + blockDim.x * blockIdx.x;
    int y = threadIdx.y + blockDim.y * blockDim.y;
    int offset = x + y * gridDim.x * blockDim.x;

    unsigned char value = src[offset];
    //原子操作，否则会出错。
    atomicAdd(&out[value],1);
}

//refer: Parallel Prefix Sum (Scan) with CUDA, Mark Harris, April 2007
//link : http: //developer.download.nvidia.com/compute/cuda/2_2/sdk/website/projects/scan/doc/scan.pdf
__global__ void scan_downsweep(unsigned int *in,unsigned int *out,int n)
{
    int Idx = threadIdx.x ;
     extern __shared__ float sdata[];
    sdata[Idx] = in[Idx];
    __syncthreads();

    unsigned int offset=1;

    for(unsigned int i=n>>1;i>0;i>>=1)  //build sum iin place up the tree
    {
        __syncthreads();
        if(Idx<i)
        {
            int ai = offset*(2*Idx+1)-1;
            int bi = offset*(2*Idx+2)-1;
            sdata[bi] += sdata[ai];
        }
        offset *= 2 ;
    }

    if(Idx==0) sdata[n-1]=0;            //travers down tree &build scan
    for(unsigned int i=1;i<n;i<<=1)
    {
        offset >>= 1;
        __syncthreads();
        if(Idx<i)
        {
            int ai = offset*(2*Idx+1)-1;
            int bi = offset*(2*Idx+2)-1;

            int tmp    =sdata[ai];
            sdata[ai]  =   sdata[bi];
            sdata[bi] +=   tmp;
        }
    }
    __syncthreads();
    //out[Idx]=sdata[Idx];
    out[Idx]=((float)sdata[Idx]/512/512*256)-1;
}
__global__ void equal_remap(char *img,char *out,unsigned int *idetity)
{
    int x = threadIdx.x+blockIdx.x*blockDim.x;
    int y = threadIdx.y+blockIdx.y*blockDim.y;
    int offset = x+y*blockDim.x*gridDim.x;

    unsigned char value = img[offset];

    out[offset] = (unsigned char) idetity[value];


}

#define DIM             512
#define GRAY_IDENTITY   256
int main(int argc,char** argv)
{

    Mat img_src = imread(IMAGE_DIR);
    Mat img_gray;
    cvtColor(img_src,img_gray,CV_RGB2GRAY);
    resize(img_gray,img_gray,Size(DIM,DIM));

    unsigned int *dev_histogram;
    cudaMalloc((void**)&dev_histogram,GRAY_IDENTITY*sizeof(int));
    cudaMemset(dev_histogram,0,GRAY_IDENTITY*sizeof(int));

    /*build the histogram of the image*/
    char *dev_gray;
    cudaMalloc((void**)&dev_gray,DIM*DIM);
    cudaMemcpy(dev_gray,img_gray.data,DIM*DIM,cudaMemcpyHostToDevice);
    dim3 blocks(DIM/16,DIM/16);dim3 threads(16,16);
    histogram_init<<<blocks,threads>>>(dev_gray,dev_histogram);
    int host_histogram[GRAY_IDENTITY]={0};
    cudaMemcpy(host_histogram,dev_histogram,GRAY_IDENTITY*sizeof(int),cudaMemcpyDeviceToHost);

    /*equalize the histogram*/
    unsigned int *dev_histogram_equal;
    cudaMalloc((void**)&dev_histogram_equal,GRAY_IDENTITY*sizeof(int));
    blocks=dim3(1,1);threads=dim3(GRAY_IDENTITY,1);
    scan_downsweep<<<blocks,threads,sizeof(int)*GRAY_IDENTITY>>>(dev_histogram,dev_histogram_equal,GRAY_IDENTITY);
    int host_histogram_equal[GRAY_IDENTITY]={0};
    cudaMemcpy(host_histogram_equal,dev_histogram_equal,GRAY_IDENTITY*sizeof(int),cudaMemcpyDeviceToHost);

    for(int i=0;i<GRAY_IDENTITY/4;i++)
    {
       //if(host_histogram[i]!=0)
         printf("bin[%3d]:%6d -> %6d| bin[%3d]:%6d -> %6d| bin[%3d]:%6d -> %6d| bin[%3d]:%6d -> %6d\n",
                4*i,host_histogram[4*i],host_histogram_equal[4*i],
                4*i+1,host_histogram[4*i+1],host_histogram_equal[4*i+1],
                4*i+2,host_histogram[4*i+2],host_histogram_equal[4*i+2],
                4*i+3,host_histogram[4*i+3],host_histogram_equal[4*i+3]);
    }

    /*remap the image use the equalized histogram*/
    char *dev_equalImg;
    cudaMalloc((void**)&dev_equalImg,DIM*DIM);
    blocks=dim3(DIM/16,DIM/16);threads=dim3(16,16);
    equal_remap<<<blocks,threads>>>(dev_gray,dev_equalImg,dev_histogram_equal);
    Mat equalImg(DIM,DIM,CV_8UC1,Scalar(255));
    cudaMemcpy(equalImg.data,dev_equalImg,DIM*DIM,cudaMemcpyDeviceToHost);




    cudaFree(dev_gray);
    cudaFree(dev_histogram);
    cudaFree(dev_histogram_equal);
    imshow("gray_img",img_gray);
    imshow("histogram equalization use GPU",equalImg);

    waitKey(0);
    return 0;
}

源码下载：GitHub
参考：Parallel Prefix Sum (Scan) with CUDA, Mark Harris, April 2007