均值滤波原理

在卷积形式的均值滤波中，需要把输入图对应模板内所有点像素点相加求平均，代替原像素点。

对于输入输出图像大小相同的均值滤波算法而言。边界的处理情况是不同的。超出图像范围的部分进行补零，这会使边界的像素值在均值滤波之后均值变小。以图一为例，左上角四个元素的均值为4。边界补零后计算3×3区域内的均值为1.8，明显均值偏小。本文算法边界采用自适应方式，比如该例子只计算四个像素的均值。

积分图

图像积分图中每个点的值是原图像中该点左上角的所有像素值之和。
$$I(x,y)=\sum _{i=0}^{i=y}\sum _{j=0}^{j=x}G(x,y)$$
I表示积分图像，G表示原始图像。

有了积分图的概念，我们可以计算图像大小范围内，任意矩形区域的像素之和。

如上图所示，要计算原图 G 中 R 区域像素之和时，可以转化为积分图I中元素的两次减法一次加法操作：$$S（R）=I4-I2-I3+I1$$如此做可以避免对原图R区域内像素累加操作，假设R区域长宽分别为m，n。计算时间复杂度可以从 O(mn) 降至 O(3)。

CPU版本一

采用卷积模板内像素值累加再求平均值的方式。

void mean_filter1(float* img_in, float*img_out, int size, int H, int W)
{
	int left = 0, right = 0, up = 0, down = 0;
	for (int i = 0; i < H; ++i)
	{
		for (int j = 0; j < W; ++j)
		{
			left = MAX(j-size, 0);
			right = MIN(j+size,W-1);
			up = MAX(i-size, 0);
			down = MIN(i+size,H-1);
			float sum = 0.0;
			for (int m = up; m <= down; m++)
			{
				for (int k = left; k <= right; k++)
				{
					sum += img_in[m*W+k];
				}
			}
			int num = (right - left + 1)*(down - up + 1);
			img_out[i*W + j] = sum / num;
		}
	}
}

CPU版本二

采用积分图的方式，先计算积分图再计算均值。

根据卷积模板大小和像素点位置确定积分图中四个顶点的位置）I1、I2、I3、I4。（A为存储积分图的矩阵，w为矩阵宽度，h为矩阵高度，x，y为待处理像素点的横纵坐标，size为模板一半长度）
$$\begin{gathered} I1=A[y-size-1][x-size-1] \\ I2=A[y-size-1][x+size] \\ I3=A[y+size][x-size-1] \\ I4=A[y+size][x+size] \end{gathered}$$
需要注意的是，本文使用的边界填充方法是超过图像边界的部分不参与均值计算，所以在边界位置需要单独处理。I1点索引如果在落在图像左或上边界之外，则I1=0；I2如果在图像上边界之外则I2=0，如果在上边界内右边界之外则I2=A[y-size-1][w-1]；I3如果在左边界之外则I3=0，如果在左边界内下边界之外则I3=A[h-1][x-size-1]；I4如果在右边界或者下边界之外，I4越界的坐标取边界值。

void mean_filter2(float* img_in, float*img_out, int size, int H, int W)
{
	/*横向相加*/
	for (int i = 0; i < H; ++i)
	{
		for (int j = 1; j < W; ++j)
		{
			img_in[i*W + j] += img_in[i*W + j - 1];
		}
	}
	/*纵向相加*/
	for (int i = 0; i < W; ++i)
	{
		for (int j = 1; j < H; ++j)
		{
			img_in[j*W + i] += img_in[(j-1)*W + i];
		}
	}

	int left = 0, right = 0, up = 0, down = 0;
	for (int i = 0; i < H; ++i)
	{
		for (int j = 0; j < W; ++j)
		{
			left = MAX((j - size-1), -1);
			right = MIN(j + size, W - 1);
			up = MAX((i - size-1), -1);
			down = MIN(i + size, H - 1);
			float point1=0, point2=0, point3 = 0;
			if (left != -1 && up != -1)
			{
				point1 = img_in[up*W + left];
				point2 = img_in[up*W + right];
				point3 = img_in[down*W + left];
			}
			if (left != -1)
			{
				point3 = img_in[down*W + left];
			}
			if (up != -1)
			{
				point2 = img_in[up*W + right];
			}
			float sum = img_in[down*W+right]
				-point2-point3+point1;
			
			int num = (right - left)*(down - up);
			img_out[i*W + j] = sum / num;
		}
	}
}

GPU版本

建立积分图时过程中，横向求和时，每一个线程计算积分图的每一行，纵向求和时，每一个线程计算图像的每一列。计算均值过程中每一个线程负责每个像素点的计算。

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include"iostream"
#include<opencv2/opencv.hpp>
#include<algorithm>
#include <windows.h>
#include <profileapi.h>
using namespace std;
using namespace cv;

__global__ void add_y(float* img_in, int W, int H)
{
	const int id = blockDim.x*blockIdx.x + threadIdx.x;
	if (id<W)
	{
		for (int j = 1; j<H; j++)
		{
			img_in[j*W + id] += img_in[(j - 1)*W + id];			
		}
	}
}

__global__ void add_x(float* img_in, int W, int H)
{
	const int id = blockDim.x*blockIdx.x + threadIdx.x;
	if (id<H)
	{
		for (int i = 1; i<W; i++)
		{
			img_in[id*W + i] += img_in[id*W + i - 1];
		}
	}
}

__global__ void mean_filter(float* img_out, float* img_in, int W, int H, int size)
{
	const int idx = blockDim.x*blockIdx.x + threadIdx.x;
	const int idy = blockDim.y*blockIdx.y + threadIdx.y;
	if (idx<W && idy<H)
	{		

		int left = MAX((idx - size - 1), -1);
		int right = MIN(idx + size, W - 1);
		int up = MAX((idy - size - 1), -1);
		int down = MIN(idy + size, H - 1);
		float point1 = 0, point2 = 0, point3 = 0;
		if (left != -1 && up != -1)
		{
			point1 = img_in[up*W + left];
			point2 = img_in[up*W + right];
			point3 = img_in[down*W + left];
		}
		else if (left != -1)
		{
			point3 = img_in[down*W + left];
		}
		else if (up != -1)
		{
			point2 = img_in[up*W + right];
		}
		float sum = img_in[down*W + right]
			- point2 - point3 + point1;

		int num = (right - left)*(down - up);
		img_out[idy*W + idx] = sum / num;
	}
}

void main()
{
	LARGE_INTEGER nFreq;
	LARGE_INTEGER nBeginTime;
	LARGE_INTEGER nEndTime;

	Mat img = imread("D:/512.jpg", CV_LOAD_IMAGE_GRAYSCALE);
	int H = img.rows;
	int W = img.cols;
	imshow("origin image", img);
	waitKey(0);

	img.convertTo(img, CV_32FC1);	
	float* img_in = NULL;	
	float* img_out = NULL;	
	cudaMalloc(&img_in, H*W * sizeof(float));
	cudaMalloc(&img_out, H*W * sizeof(float));
	cudaMemcpy(img_in, img.data, H*W * sizeof(float), cudaMemcpyHostToDevice);

	QueryPerformanceFrequency(&nFreq);
	QueryPerformanceCounter(&nBeginTime);			
	//纵向相加
	dim3 blocksize(16, 1);
	dim3 gridsize((W + blocksize.x - 1) / blocksize.x, 1);
	add_y << <gridsize, blocksize >> >(img_in, W, H);
	//横向相加
	blocksize = dim3(16, 1);
	gridsize = dim3((H + blocksize.x - 1) / blocksize.x, 1);
	add_x << <gridsize, blocksize >> >(img_in, W, H);
	//滤波
	blocksize = dim3(16, 16);
	gridsize = dim3((W + blocksize.x - 1) / blocksize.x,
		(H + blocksize.y - 1) / blocksize.y);
	mean_filter << <gridsize, blocksize >> >(img_out, img_in, W, H, 2);
	
	QueryPerformanceCounter(&nEndTime);
	double time = (double)(nEndTime.QuadPart - nBeginTime.QuadPart) / (double)nFreq.QuadPart;
	cout << time << "s" << endl;

	cudaMemcpy(img.data, img_out, H*W * sizeof(float), cudaMemcpyDeviceToHost);
	cudaFree(img_in);
	cudaFree(img_out);
	

	img.convertTo(img, CV_8UC1);
	imshow("blur image", img);
	waitKey(0);
}

版本比较

本次实验代码运行环境为CPU：intel i5-9300H，GPU：NVIDIA GTX1660ti，测试时间去除了内存分配占用的时间，仅统计原图矩阵到均值滤波后的矩阵过程的时间。

算法分析

CPU1方案：假设图像为边长为n的方阵，卷积模板为边长为m的正方形区域，对于n2个像素点，每一个像素点周围m2个像素进行累加操作，算法的时间复杂度为O（n2m2）。从运行结果看，当图像的大小加倍，卷积核大小不变时时，计算时间翻了四倍。当图像比例不变卷积核大小翻2.2倍时，计算时间翻了4.4倍左右。当图像大小翻倍和卷积核大小为2.2倍时，计算时间翻了19倍左右。这与理论上算法的时间复杂度一致。
CPU2方案：建立积分图时，每一行进行n-1次相加，一共n行，每一列也是如此，时间复杂度O(n2)。滤波时，对于每一个像素点执行3次加减操作，一共有n2个像素点，时间复杂度为O(n2)。从运算结果看，图像大小加倍时，计算时间翻了4倍左右，卷积核的变化对该方案没有影响，与理论一致。
GPU方案：建表过程，每一个线程进行n-1次加法操作。假设每个线程同时运行，计算时间为O(n)，滤波过程中，对于每一个线程，只用处理一个像素点，每个像素点进行一次3次加减操作计算时间为O(1)。总体时间复杂度为O(n)。按理论讲图像大小加倍，计算时间应该加倍，可能由于时间太短，测量精度没有达到。