作者：翟天保Steven
版权声明：著作权归作者所有，商业转载请联系作者获得授权，非商业转载请注明出处
实现原理
基于累计直方图的中值滤波算法详情参见：
OpenCV-基于累计直方图的中值滤波算法-CSDN博客
为了进一步提升算法性能，我尝试采用CUDA进行提速，但在实际实现中发现，如果在CUDA的核函数中频繁申请释放内存，其速度还不如CPU。因此再结合共享内存实现，在GPU编程中，共享内存是每个线程块共享的内存空间，它的读写速度比全局内存快得多。
功能函数代码
// 在GPU上执行中值滤波
__global__ void filterMedian_CUDA(const uchar* input, uchar* output, int rows, int cols, int filterWindowSize)
{
__shared__ int sharedHis[256 * TILE_WIDTH * TILE_WIDTH];
int i = blockIdx.y * blockDim.y + threadIdx.y;
int j = blockIdx.x * blockDim.x + threadIdx.x;
if (i < rows && j < cols)
{
int r = filterWindowSize / 2;
int ms = max(i - r, 0);
int me = min(i + r, rows - 1);
int ns = max(j - r, 0);
int ne = min(j + r, cols - 1);
// 初始化
int count = 0;
int idx = ((threadIdx.y * blockDim.x + threadIdx.x)) * 256;;
for (int k = 0; k < 256; ++k)
{
sharedHis[idx + k] = 0;
}
// 直方图累计求中值
for (int m = ms; m <= me; ++m)
{
for (int n = ns; n <= ne; ++n)
{
sharedHis[idx + input[m * cols + n]]++;
count++;
}
}
int count_ = 0;
int thresh = count / 2 + 1;
int mid = 0;
for (int k = 0; k < 256; ++k)
{
count_ += sharedHis[idx + k];
if (count_ >= thresh)
{
mid = k;
break;
}
}
output[i * cols + j] = mid;
}
}C++测试代码
Test.h
#include <math.h>
#include <vector>
#include <iostream>
#include <algorithm>
#include <opencv2/opencv.hpp>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
using namespace std;
using namespace cv;
#define TILE_WIDTH 4
// 预准备过程
void warmupCUDA();
// 中值滤波-CUDA
cv::Mat filterMedian_CUDA(cv::Mat input, int filterWindowSize);
#endif // TEST_HTest.cu
#include "Test.h"
// 预准备过程
void warmupCUDA()
{
float *dummy_data;
cudaMalloc((void**)&dummy_data, sizeof(float));
cudaFree(dummy_data);
}
// 在GPU上执行中值滤波
__global__ void filterMedian_CUDA(const uchar* input, uchar* output, int rows, int cols, int filterWindowSize)
{
__shared__ int sharedHis[256 * TILE_WIDTH * TILE_WIDTH];
int i = blockIdx.y * blockDim.y + threadIdx.y;
int j = blockIdx.x * blockDim.x + threadIdx.x;
if (i < rows && j < cols)
{
int r = filterWindowSize / 2;
int ms = max(i - r, 0);
int me = min(i + r, rows - 1);
int ns = max(j - r, 0);
int ne = min(j + r, cols - 1);
// 初始化
int count = 0;
int idx = ((threadIdx.y * blockDim.x + threadIdx.x)) * 256;;
for (int k = 0; k < 256; ++k)
{
sharedHis[idx + k] = 0;
}
// 直方图累计求中值
for (int m = ms; m <= me; ++m)
{
for (int n = ns; n <= ne; ++n)
{
sharedHis[idx + input[m * cols + n]]++;
count++;
}
}
int count_ = 0;
int thresh = count / 2 + 1;
int mid = 0;
for (int k = 0; k < 256; ++k)
{
count_ += sharedHis[idx + k];
if (count_ >= thresh)
{
mid = k;
break;
}
}
output[i * cols + j] = mid;
}
}
// 中值滤波-CUDA
cv::Mat filterMedian_CUDA(cv::Mat input, int filterWindowSize)
{
int row = input.rows;
int col = input.cols;
int size = row * col * sizeof(uchar);
uchar *d_input, *d_output;
// 在GPU上分配内存
cudaMalloc((void**)&d_input, size);
cudaMalloc((void**)&d_output, size);
// 将输入数据从主机内存复制到GPU内存
cudaMemcpy(d_input, input.data, size, cudaMemcpyHostToDevice);
// 计算CUDA核函数执行的网格和块大小
dim3 threadsPerBlock(TILE_WIDTH, TILE_WIDTH);
dim3 blocksPerGrid((col + threadsPerBlock.x - 1) / threadsPerBlock.x, (row + threadsPerBlock.y - 1) / threadsPerBlock.y);
// 在GPU上执行中值滤波
filterMedian_CUDA << <blocksPerGrid, threadsPerBlock >> > (d_input, d_output, row, col, filterWindowSize);
// 将结果从GPU内存复制回主机内存
cv::Mat output(row, col, CV_8UC1);
cudaMemcpy(output.data, d_output, size, cudaMemcpyDeviceToHost);
// 释放GPU内存
cudaFree(d_input);
cudaFree(d_output);
return output;
}
main.cpp
#include <time.h>
#include <omp.h>
#include <vector>
#include <iostream>
#include <algorithm>
#include <opencv2/opencv.hpp>
#include "Test.h"
using namespace std;
// 基于累计直方图的中值滤波算法
cv::Mat filterMedian_CPU(cv::Mat input, int filterWindowSize)
{
int row = input.rows;
int col = input.cols;
int r = filterWindowSize / 2;
cv::Mat output(row, col, CV_8UC1);
#pragma omp parallel for
for (int i = 0; i < row; ++i)
{
for (int j = 0; j < col; ++j)
{
int ms = std::max(i - r, 0);
int me = std::min(i + r, row - 1);
int ns = std::max(j - r, 0);
int ne = std::min(j + r, col - 1);
// 初始化
int count = 0;
int histogram[256] = { 0 };
// 直方图累计求中值
for (int m = ms; m <= me; ++m)
{
for (int n = ns; n <= ne; ++n)
{
histogram[input.at<uchar>(m, n)]++;
count++;
}
}
int count_ = 0;
int thresh = count / 2 + 1;
int mid = 0;
for (int k = 0; k < 256; ++k)
{
count_ += histogram[k];
if (count_ >= thresh)
{
mid = k;
break;
}
}
output.at<uchar>(i, j) = mid;
}
}
return output;
}
int main(void)
{
// 预准备过程
warmupCUDA();
// 测试
cv::Mat input = cv::imread("test.jpg", 0);
cv::Mat output0, output1;
int filterWindowSize = 15;
cout << "filterWindowSize:" << filterWindowSize << endl;
cout << "size: " << input.cols << " * " << input.rows << endl;
// 在CPU上执行中值滤波-累计直方图法
clock_t sc, ec;
sc = clock();
output0 = filterMedian_CPU(input, filterWindowSize);
ec = clock();
cout << "His CPU time:" << double(ec - sc) / 1000 << endl;
// 在GPU上执行中值滤波-累计直方图法
clock_t sg, eg;
sg = clock();
output1 = filterMedian_CUDA(input, filterWindowSize);
eg = clock();
cout << "His GPU time:" << double(eg - sg) / 1000 << endl;
// 检查结果
int row = input.rows;
int col = input.cols;
bool flag = true;
for (int i = 0; i < row; ++i)
{
for (int j = 0; j < col; ++j)
{
if (output0.at<uchar>(i, j) != output1.at<uchar>(i, j))
{
flag = false;
break;
}
}
if (!flag)
{
break;
}
}
if (flag)
{
std::cout << "ok!" << std::endl;
}
else
{
std::cout << "ng!" << std::endl;
}
cv::Mat test0 = output0.clone();
cv::Mat test1 = output1.clone();
cout << "CUDA test." << endl;
return 0;
}测试效果
如上图所示，分别是原图、CPU结果和GPU结果，在速度方面，对1920*1080的图像，在窗口尺寸为15*15时，我的电脑运行速度分别是0.2s和0.1s。
我在不同窗尺寸下进行了几轮测试，测试结果仅供参考。
在CUDA编程中，内存访问是一个关键的性能瓶颈。频繁地申请和释放内存会增加内存管理的开销，这可能会导致程序性能下降。为了提高CUDA程序的效率，可以尽量减少内存分配和释放的次数。一种常见的方法是预先分配一定大小的内存，并在程序执行过程中重复使用这些内存。此外，可以考虑使用共享内存等技术来减少对全局内存的访问，以提高访存效率。
如果函数有什么可以改进完善的地方，非常欢迎大家指出，一同进步何乐而不为呢~
如果文章帮助到你了，可以点个赞让我知道，我会很快乐~加油！
