使用sycl实现对图像高斯模糊处理

WantCc

已于 2023-11-23 11:13:07 修改

阅读量122

点赞数 1

分类专栏： sycl 文章标签： c++ 图像处理

于 2023-11-23 11:12:14 首次发布

本文链接：https://blog.csdn.net/qq_36912307/article/details/134477790

版权

sycl 专栏收录该内容

2 篇文章 0 订阅

订阅专栏

一、sycl

SYCL是一个免版税的跨平台抽象层，基于OpenCL的基本概念，可移植性和效率，使得异构处理器的代码可以使用完全标准的“单一来源”风格编写C ++。 SYCL支持单一源代码开发，其中C ++模板函数可以包含主机代码和设备代码，以构建使用OpenCL加速的复杂算法，然后在不同类型数据的源代码中重复使用它们。

二、高斯模糊

高斯模糊，也叫高斯平滑，是在Adobe Photoshop、GIMP以及Paint.NET等图像处理软件中广泛使用的处理效果，通常用它来减少图像噪声以及降低细节层次。这种模糊技术生成的图像，其视觉效果就像是经过一个半透明屏幕在观察图像，这与镜头焦外成像效果散景以及普通照明阴影中的效果都明显不同。最常见的就是当你把手机状态栏下拉的时候可以看到桌面的画面作为背景，这个时候的左面就是被高斯模糊的。这样可以给人机交互体验上有一个很大的提升。

从数学的角度来看，图像的高斯模糊过程就是图像与正态分布做卷积。由于正态分布又叫作“高斯分布”，所以这项技术就叫作高斯模糊。图像与圆形方框模糊做卷积将会生成更加精确的焦外成像效果。由于高斯函数的傅立叶变换是另外一个高斯函数，所以高斯模糊对于图像来说就是一个低通滤波器。在二维空间定义为：

其中r是模糊半径，σ是正态分布的标准偏差。

分布不为零的像素组成的卷积矩阵与原始图像做变换。每个像素的值都是周围相邻像素值的加权平均。原始像素的值有最大的高斯分布值，所以有最大的权重，相邻像素随着距离原始像素越来越远，其权重也越来越小。这样进行模糊处理比其它的均衡模糊滤波器更高地保留了边缘效果。简单来说就是可以通过上面那个公式来计算出来高斯核，当参数确定之后，高斯核也是确定的。然后我们就可以将目标图像与计算好的正态分布做卷积，已选好的目标像素值为中心，根据高斯卷积核对周围的像素值进行加权求和。下面就是一个计算σ = 0.84089642的高斯分布生成的示例矩阵。注意中心元素 (4,4)处有最大值，随着距离中心越远数值对称地减小。

旁边是卷积操作动态图，与高斯模糊实现的过程类似，可以借助这个动画来理解高斯模糊的过程

四、代码实现及相关sycl知识讲解

本来想单独开一篇文章来讲解sycl使用的，但是博主对sycl的了解也不是很多，大部分都是基于这次使用才了解到的，所以就简单一块讲了。

首先需要是将图像载入，这里用到的是stb_image库GitHub - nothings/stb: stb single-file public domain libraries for C/C++；头文件以及图像载入实现如下

#define STB_IMAGE_IMPLEMENTATION
#include "stb/stb_image.h"
#define STB_IMAGE_WRITE_IMPLEMENTATION
#include "stb/stb_image_write.h"
 int inputWidth, inputHeight, inputChannels;
  const int numChannels = 4;
  void* inputData = nullptr;
  void* outputData = nullptr；
  inputData = stbi_load(argv[1], &inputWidth, &inputHeight, &inputChannels,numChannels);//加载图像
  outputData = new char[inputWidth * inputHeight * numChannels];//设置输出图像

数据并行性--一些数据结构概念

在 SYCL中， range对象用于描述一维空间的范围。其中range 类表示一个范围，包含一个整数值，用于表示该范围的大小。它可以在 SYCL 内核函数中与 nd_range 一起使用，指定并行计算的执行范围。SYCL 还提供了 range 的二维和三维版本，分别是 range<2> 和 range<3>，用于描述二维和三维空间的范围。需要注意的是range 主要用于表示全局的大小，而 nd_range 则是一个更复杂的对象，除了全局大小外，还包含局部大小，用于更精确地描述并行计算的执行范围和分配。nd_range 可以在 SYCL 内核函数中与 item 或 nd_item 一起使用，以确定每个工作项的位置和范围。

image<2> 是 SYCL中的一个类，用于表示二维图像数据。image<2> 类提供了图像数据的访问和操作接口，可以在 SYCL 内核函数中使用。

inputData 是一个用于存储图像数据的缓冲区或指针。
co::rgba 是一个像素格式，表示图像中每个像素的通道顺序为 RGBA。
ct::unorm_int8 是一个像素数据类型（Channel Type），表示每个通道使用 8 位无符号整数。
imgRange 是一个 range<2> 对象，用于指定图像的宽度和高度。

  range<2> imgRange(inputHeight, inputWidth);
  range<2> gaussianRange( guasize*stddev+1,  guasize*stddev+1);
 image<2> image_in(inputData, co::rgba, ct::unorm_int8, imgRange);
 image<2> image_out(outputData, co::rgba, ct::unorm_int8, imgRange);

前面提到nd_range 则是一个更复杂的对象，除了全局大小外，还包含局部大小。其对应的构造参数如下

nd_range(range<dimensions> globalSize, range<dimensions> localSize);

globalSize：全局范围的大小，表示并行计算的总体范围。这通常是一个 range 对象，可以是一维、二维或三维的。
localSize：局部范围的大小，表示并行计算中每个工作组的大小。这通常是一个 range 对象，可以是一维、二维或三维的。如果不需要使用局部范围，则可以将其设置为默认构造函数创建的空范围。

queue

简单理解，就是我们创建一个任务队列，然后将我们的任务提交到任务队列，然后执行

queue Q;
Q.submit(
  /* device code action */
);

我们也可以选择特定类别的设备创建队列来执行任务，一个queue映射到一台device：映射发生在队列对象构造时，并且随后无法更改。

选择特定类别的设备
queue Q_cpu{default_selector{}};
queue Q_cpu{cpu_selector{}};
queue Q_device{gpu_selector{}};
queue Q_accelerator{accelerator_selector{}};
auto Q = queue{my_selector{}};

所提供的handler（handler &cgh）可以更好地帮助我们控制代码如何提交到队列。

auto Q = queue{my_selector{}};

Q.submit([&](handler &cgh){
主机代码，用于设置任务图中相应节点的依赖关系。 主机代码在提交后立即执行。 

 该操作在设备上异步执行。 此外，并行工作操作还需要一个执行范围和一个内核函数。
 cgh.parallel_for(range<1>{sz}, [=](auto &idx){
    /* kernel code */
 });
});

在一般情况下，使用 CPU 上的时钟观测 SYCL 内核函数的运行时间是不准确的。这是因为 SYCL 内核函数通常在加速器设备上执行，而不是在主机的 CPU 上执行。这与opencl中思路一样

我们在openCL中为了查看在加速设备中处理的时间需要这样操作

（有需要的老哥自取）

clCreateCommandQueue第三个参数一定要设置上ENABLE，都是泪，然后就可以使用Get获取对应的时间戳了

cl_command_queue queue= clCreateCommandQueue(context, device, CL_QUEUE_PROFILING_ENABLE, NULL);
clWaitForEvents(1, &event);
cl_ulong start, end;
clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &start, nullptr);
clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_END, sizeof(cl_ulong), &end, nullptr);
unsigned long elapsed = (unsigned long)(end - start);
printf("\nExecution time in milliseconds = %.3f ms\n", elapsed/1000000.0  );

在sycl中也有类似的用法，submit函数返回一个sycl::event 对象，它代表了任务的执行事件。可以使用该事件对象来查询任务的执行状态、等待任务完成或设置事件之间的依赖关系。然后调用队列的 wait_and_throw 函数来等待之前提交的任务执行完成。wait_and_throw 函数会阻塞当前线程，直到队列中的所有任务都执行完成。

static void ReportTime(const std::string &msg, sycl::event e) {
  cl_ulong time_start =
      e.get_profiling_info<sycl::info::event_profiling::command_start>();

  cl_ulong time_end =
      e.get_profiling_info<sycl::info::event_profiling::command_end>();

  double elapsed = (time_end - time_start) / 1e6;
  std::cout << msg << elapsed << " milliseconds\n";
}
sycl::event e1 = myQueue.submit([&](handler& cgh) {
/*
            要执行的操作
*/
}
myQueue.wait_and_throw();
ReportTime("时间为：",e1);

使用buffer和accessor进行数据管理

buffer是 1、2 或 3 维内存位置的句柄。它指定内存位置以及可以访问的位置：host、device或两者。因此，buffer不拥有内存：它只是内存的一个受限视图。我们不直接处理buffer，而是使用accessor，accessor是用于在内核函数中对buffer进行访问和操作的对象。accessor提供了一种机制，使得内核函数可以读取和写入buffer中的数据。sycl还提供了images，images与buffer类型类似，同时为图像处理量身定制的额外功能。

buffer“告诉”运行时数据如何布局，而accessor“告诉”它我们将如何读取和写入底层内存。

accessor可以通过多个参数进行模板化，（有需要的同学可以去看一下这个文档https://www.khronos.org/files/sycl/sycl-2020-reference-guide.pdf）不过我们这里只掌握一种最简单的使用方法就可以

即先创建一个buffer，然后传入sccessor 设置读写模式，例如下面：

auto inBuf = sycl::buffer{inputImage.data(), inBufRange};
auto outBuf = sycl::buffer<float, 2>{outBufRange};
auto filterBuf = sycl::buffer{filter.data(), filterRange};
sycl::accessor inputAcc{inBuf, cgh, sycl::read_only};
sycl::accessor outputAcc{outBuf, cgh, sycl::write_only};
sycl::accessor filterAcc{filterBuf, cgh, sycl::read_only};

image 类用于表示图像对象，而 sampler 类用于表示采样器对象。通过将图像对象和采样器对象传递给内核函数，可以在内核中使用sampler 来对图像进行采样操作。

在 SYCL 中，sampler 类是用于指定图像采样器对象，在内核中对图像进行采样操作。它与 image 类一起使用，指定在图像中进行采样时的采样参数，例如过滤器类型、边界处理等。

使用方法如下，具体参数设置可见DPC++ Runtime: include/sycl/sampler.hpp Source File

sycl::sampler smp(coordinate_normalization_mode::unnormalized, addressing_mode::clamp, filtering_mode::nearest);

基本工作做的差不多了，就可以将任务提交到device上进行异步执行了。

#include <CL/sycl.hpp>
#include <vector>
#define _USE_MATH_DEFINES 
#include <cmath>
#include <iostream>
#ifdef _MSC_VER
typedef unsigned int uint;
#endif
#define STB_IMAGE_IMPLEMENTATION
#include "stb/stb_image.h"
#define STB_IMAGE_WRITE_IMPLEMENTATION
#include "stb/stb_image_write.h"
#include "dpc_common.hpp"
class fillGaussian;
class GaussianKernel;
      float sum ;

using namespace cl::sycl;
using co = cl::sycl::image_channel_order;
using ct = cl::sycl::image_channel_type;
void gaussianBlur(unsigned char* image, int width, int height, int channels, float sigma);
range<2> get_optimal_local_range(cl::sycl::range<2> globalSize,cl::sycl::device d) {
    range<2> optimalLocalSize(0,0); 
  if (d.is_gpu()) {
    optimalLocalSize = range<2>(64, 1);
      std::cout<<"this is GPU"<<std::endl;
     
  } else {
    optimalLocalSize = range<2>(64, 1);
  }
  for (int i = 0; i < 2; ++i) {
    while (globalSize[i] % optimalLocalSize[i]) {
      optimalLocalSize[i] = optimalLocalSize[i] >> 1;
    }
  }
  return optimalLocalSize;
}
static void ReportTime(const std::string &msg, sycl::event e) {
  cl_ulong time_start =
      e.get_profiling_info<sycl::info::event_profiling::command_start>();

  cl_ulong time_end =
      e.get_profiling_info<sycl::info::event_profiling::command_end>();

  double elapsed = (time_end - time_start) / 1e6;
  std::cout << msg << elapsed << " milliseconds\n";
}
int main(int argc, char* argv[]) {

  int inputWidth, inputHeight, inputChannels;
  const int numChannels = 4;
  void* inputData = nullptr;
  void* outputData = nullptr;
  if (argc < 2) {
    std::cout
        << "没有图片"
        << std::endl;
  }
  inputData = stbi_load(argv[1], &inputWidth, &inputHeight, &inputChannels,numChannels);//加载图像  
  if (inputData == nullptr) {
    std::cout << "Failed to load image file (is argv[1] a valid image file?)"<< std::endl;
    return 1;
  }
  outputData = new char[inputWidth * inputHeight * numChannels];//设置输出图像
  const float pi = M_PI;
  static constexpr auto stddev = 2;
  const int guasize = 6;
  range<2> imgRange(inputHeight, inputWidth);
  range<2> gaussianRange( guasize*stddev+1,  guasize*stddev+1);
  float(*printGuaRan)[guasize*stddev+1] = new float[guasize*stddev+1][guasize*stddev+1];
     auto prop_list = property_list{property::queue::enable_profiling()};
     default_selector d_selector;
    auto platforms =  sycl::platform::get_platforms();
    queue myQueue(d_selector,dpc_common::exception_handler, prop_list);
    std::cout << "Device: " << myQueue.get_device().get_info<info::device::name>() << "\n";
  {
   buffer gaussian(reinterpret_cast<float *>(printGuaRan), range(guasize*stddev+1, guasize*stddev+1));
        for(int i=0;i<guasize*stddev+1;i++){
        for (int j=0;j<guasize*stddev+1;j++){
             printGuaRan[i][j] = (1/(2 * M_PI ))*std::exp(-1.f * (i * i + j * j) / 2) ;
            sum+=printGuaRan[i][j];
        }
    }
    image<2> image_in(inputData, co::rgba, ct::unorm_int8, imgRange);
    image<2> image_out(outputData, co::rgba, ct::unorm_int8, imgRange);
    sycl::event e1 = myQueue.submit([&](handler& cgh) {
      auto r = get_optimal_local_range(imgRange, myQueue.get_device());
      auto myRange = nd_range<2>(imgRange, r);
      accessor<float4, 2, access::mode::read, access::target::image> inPtr(image_in, cgh);
      accessor<float4, 2, access::mode::write, access::target::image> outPtr(image_out, cgh);
      auto globalGaussian = gaussian.get_access<access::mode::read>(cgh);
      sampler smpl(coordinate_normalization_mode::unnormalized,addressing_mode::clamp_to_edge,filtering_mode::nearest);
      cgh.parallel_for<GaussianKernel>(myRange, [=](nd_item<2> itemID) {
        float4 newPixel = float4(0.0f, 0.0f, 0.0f, 0.0f);
        constexpr auto offset = guasize * stddev/2;
        for (int x = -offset; x < offset; x++) {
          for (int y = -offset; y < offset; y++) {
            auto inputCoords =
                int2(itemID.get_global_id(0)+y , itemID.get_global_id(1)+x );
            newPixel += inPtr.read(inputCoords, smpl)*(globalGaussian[y + offset][x + offset]/0.48926);
          }
        }
        auto outputCoords =
            int2(itemID.get_global_id(0), itemID.get_global_id(1));
        newPixel.w() = 1.f;
        outPtr.write(outputCoords, newPixel);
      });
    });
    myQueue.wait_and_throw();
      std::cout<<"图像的大小为:"<<inputWidth<<"*"<<inputHeight<<"*"<<inputChannels<<std::endl;
        ReportTime("时间为：",e1);
  }
    std::cout<<"********串行设计*****************";
        
    int width, height, channels;
    unsigned char* image = stbi_load(argv[1], &width, &height, &channels, 0);
    if (image == nullptr)
    {
        std::cout << "Failed to load image." << std::endl;
        return 1;
    }
    float sigma = 1.5f; // 高斯核的标准差
    
    gaussianBlur(image, width, height, channels, sigma);
    stbi_write_jpg("./output.jpg", width, height, channels, image, 100);
    stbi_image_free(image);  

  std::string outputFilePath;
  std::string inputName(argv[1]);
  auto pos = inputName.find_last_of(".");
  if (pos == std::string::npos) {
    outputFilePath = inputName + "-blurred";
  } else {
    std::string ext = inputName.substr(pos, inputName.size() - pos);
    inputName.erase(pos, inputName.size());
    outputFilePath = inputName + "-blurred" + ext;
  }
  stbi_write_png("/home/u205092/HPC/gaussian/gaussian.jpg", inputWidth, inputHeight, numChannels,
                 outputData, 0);

  std::cout << "Image successfully blurred!\n";
  return 0;
}
void gaussianBlur(unsigned char* image, int width, int height, int channels, float sigma)
{
    unsigned char* blurredImage = new unsigned char[width * height * channels];
    int kernelSize = 2 * ceil(2 * sigma) + 1;
    float* kernel = new float[kernelSize * kernelSize];
    float kernelSum = 0.0f;
    int kernelRadius = kernelSize / 2;
    for (int i = -kernelRadius; i <= kernelRadius; ++i)
    {
        for (int j = -kernelRadius; j <= kernelRadius; ++j)
        {
            float weight = exp(-(i * i + j * j) / (2 * sigma * sigma));
            kernel[(i + kernelRadius) * kernelSize + (j + kernelRadius)] = weight;
            kernelSum += weight;
        }
    }
    for (int i = 0; i < kernelSize; ++i)
    {
        for (int j = 0; j < kernelSize; ++j)
        {
            kernel[i * kernelSize + j] /= kernelSum;
        }
    }
         time_t S_begin_t  = clock();

    for (int y = 0; y < height; ++y)
    {
        for (int x = 0; x < width; ++x)
        {
            for (int c = 0; c < channels; ++c)
            {
                float blurredPixel = 0.0f;
                for (int i = -kernelRadius; i <= kernelRadius; ++i)
                {
                    for (int j = -kernelRadius; j <= kernelRadius; ++j)
                    {
                        int neighborX = std::clamp(x + j, 0, width - 1);
                        int neighborY = std::clamp(y + i, 0, height - 1);
                        blurredPixel += kernel[(i + kernelRadius) * kernelSize + (j + kernelRadius)] *
                                        image[(neighborY * width + neighborX) * channels + c];
                    }
                }
                blurredImage[(y * width + x) * channels + c] = static_cast<unsigned char>(blurredPixel);
            }
        }
    }
    time_t S_end_t  = clock();
    std::cout<<"串行计算耗费时间为： " << (double )(S_end_t - S_begin_t )/CLOCKS_PER_SEC*1000 <<" ms"<<std::endl;
    std::memcpy(image, blurredImage, width * height * channels);
    delete[] blurredImage;
    delete[] kernel;
}

编译（局限在devcloud，感觉可以将这个这个目录下的文件转移到自己的linux设备，查看不同的设备下加速效果，不过时间有限，然后高斯模糊需要的算力也不是很大，intel的这个devcloud给到的设备已经够用了，就没尝试）

icpx -I/opt/intel/oneapi/compiler/latest/linux/include/sycl -fsycl source.cpp

查看运行效果

这里设置串行处理用于计算加速比，查看加速效果。在代码中 gaussianBlur（）方法中实现。

WantCc

关注

1
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
使用sycl实现对图像高斯模糊处理

最常见的就是当你把手机状态栏下拉的时候可以看到桌面的画面作为背景，这个时候的左面就是被高斯模糊的。每个像素的值都是周围相邻像素值的加权平均。SYCL是一个免版税的跨平台抽象层，基于OpenCL的基本概念，可移植性和效率，使得异构处理器的代码可以使用完全标准的“单一来源”风格编写C ++。（局限在devcloud，感觉可以将这个这个目录下的文件转移到自己的linux设备，查看不同的设备下加速效果，不过时间有限，然后高斯模糊需要的算力也不是很大，intel的这个devcloud给到的设备已经够用了，就没尝试）
复制链接

扫一扫