拉普拉斯算子(Laplacian)可应用到图像边缘检测中。在OpenCV中当kernel大小为3*3时，支持两种kernel算子，分别为：

在OpenCV中默认的计算方式如下，假设有一个5*5的小图像，原始值依次为1,2,…25,如下图红色部分，首先将5*5映射到(5+3-1)*(5+3-1)大小，然后和3*3的kernel做累积和，因为计算结果有可能超出有效值[0, 255]范围，因此在最终还需要判断使其在有效范围内，即小于0为0，大于255为255：

如坐标为(0,0)的计算过程为：12 = 7*0+6*1+7*0+2*1+1*(-4)+2*1+7*0+6*1+7*0.

以下code分别采用两种方式实现，一种是从OpenCV中提取的code(Laplacian_函数)，一种是按照上面说明的方式实现的(Laplacian函数)，两种方式实现的结果完全一致，但第二种方式更容易理解：

laplacian.cpp:#include "funset.hpp"

#include

#include

#include

#include

#include

#include "common.hpp"

namespace {

typedef struct Rect {

int x, y, width, height;

} Rect;

typedef struct Size {

int width, height;

} Size;

typedef struct Point {

int x, y;

} Point;

int borderInterpolate(int p, int len)

{

if ((unsigned)p

;

} else {

if (len == 1) return 0;

int delta = 1;

do {

if (p

else p = len - 1 - (p - len) - delta;

} while ((unsigned)p >= (unsigned)len);

}

return p;

}

inline unsigned char saturate_cast(float v)

{

int iv = (int)(v + (v >= 0 ? 0.5f : -0.5f));

return (unsigned char)((unsigned)v <= UCHAR_MAX ? v : v > 0 ? UCHAR_MAX : 0);

}

void filter2D(const unsigned char** src, unsigned char* dst, int dststep, int count, int width, int ksize)

{

std::vector coords;

std::vector coeffs;

if (ksize == 1) {

coords = { { 1, 0 }, { 0, 1 }, { 1, 1 }, { 2, 1 }, { 1, 2 } }; // kernel non zero position: (x, y)

coeffs = { 1.f, 1.f, -4.f, 1.f, 1.f }; // kernel non zero value: 1, 1, -4, 1, 1

} else {

coords = { { 0, 0 }, { 2, 0 }, { 1, 1 }, { 0, 2 }, { 2, 2 } }; // kernel non zero position: (x, y)

coeffs = { 2.f, 2.f, -8.f, 2.f, 2.f }; // kernel non zero value: 2, 2, -8, 2, 2

}

std::vector ptrs(coords.size());

float _delta{ 0.f };

const Point* pt = &coords[0];

const float* kf = (const float*)&coeffs[0];

const unsigned char** kp = (const unsigned char**)&ptrs[0];

int nz = (int)coords.size();

for (; count > 0; count--, dst += dststep, src++) {

unsigned char* D = (unsigned char*)dst;

for (int k = 0; k

kp[k] = (const unsigned char*)src[pt[k].y] + pt[k].x;

for (int i = 0; i

float s0 = _delta;

for (int k = 0; k

s0 += kf[k] * kp[k][i];

D[i] = saturate_cast(s0);

}

}

}

int Laplacian_(const unsigned char* src_, unsigned char* dst_, int width_, int height_, int ksize_)

{

const unsigned char* src = src_;

unsigned char* dst = dst_;

const Size ksize{ 3, 3 };

const int maxBufRows = ksize.height + 3;

const Point anchor{ 1, 1 };

const Rect roi{ 0, 0, width_, height_ };

const int dx1{ 1 }, dx2{ 1 };

int borderLength = std::max(ksize.width - 1, 1);

std::vector borderTab(borderLength);

borderTab[0] = borderInterpolate(-dx1, width_);

borderTab[1] = borderInterpolate(width_, width_);

std::vector rows(maxBufRows);

const int* btab = &borderTab[0];

int srcstep{ width_ }, dststep{ width_ };

std::vector ringBuf((width_ + ksize.width - 1) * maxBufRows, 0);

int bufStep{ width_ + ksize.width - 1 };

int startY = std::max(roi.y - anchor.y, 0), startY0 = startY, rowCount{ 0 }, dstY{ 0 };

int endY = std::min(roi.y + roi.height + ksize.height - anchor.y - 1, height_);

int esz = 1;

unsigned char** brows = &rows[0];

int bufRows = (int)rows.size();

int kwidth = ksize.width;

int kheight = ksize.height, ay = anchor.y;

int _dx1 = dx1, _dx2 = dx2;

int width1 = roi.width + kwidth - 1;

int dy = 0, i = 0;

int count = endY - startY;

for (;; dst += dststep * i, dy += i) {

int dcount = bufRows - ay - startY - rowCount + roi.y;

dcount = dcount > 0 ? dcount : bufRows - kheight + 1;

dcount = std::min(dcount, count);

count -= dcount;

for (; dcount-- > 0; src += srcstep) {

int bi = (startY - startY0 + rowCount) % bufRows;

unsigned char* brow = &ringBuf[0] + bi*bufStep;

unsigned char* row = brow;

if (++rowCount > bufRows) {

--rowCount;

++startY;

}

memcpy(row + _dx1*esz, src, (width1 - _dx2 - _dx1)*esz);

for (i = 0; i

row[i] = src[btab[i]];

for (i = 0; i

row[i + (width1 - _dx2)*esz] = src[btab[i + _dx1*esz]];

}

int max_i = std::min(bufRows, roi.height - (dstY + dy) + (kheight - 1));

for (i = 0; i

int srcY = borderInterpolate(dstY + dy + i + roi.y - ay, height_);

if (srcY

if (srcY >= startY + rowCount) break;

int bi = (srcY - startY0) % bufRows;

brows[i] = &ringBuf[0] + bi*bufStep;

}

if (i

i -= kheight - 1;

filter2D((const unsigned char**)brows, dst, dststep, i, roi.width, ksize_);

}

dstY += dy;

if (dstY > roi.height) return -1;

return 0;

}

int Laplacian(const unsigned char* src_, unsigned char* dst_, int width_, int height_, int ksize_)

{

const int kernel_size{ 3 };

std::vector kernel;

if (ksize_ == 1) kernel = { 0.f, 1.f, 0.f, 1.f, -4.f, 1.f, 0.f, 1.f, 0.f };

else kernel = { 2.f, 0.f, 2.f, 0.f, -8.f, 0.f, 2.f, 0.f, 2.f };

int new_width = width_ + kernel_size - 1, new_height = height_ + kernel_size - 1;

std::unique_ptr data(new unsigned char[new_width * new_height]);

unsigned char* p = data.get();

for (int y = 0; y

if (y != 0 && y != new_height - 1) {

for (int x = 0; x

if (x == 0) {

p[y * new_width + x] = src_[(y - 1) * width_ + 1];

} else if (x == new_width - 1) {

p[y * new_width + x] = src_[(y - 1) * width_ + (width_ - 1 - 1)];

} else {

p[y * new_width + x] = src_[(y - 1) * width_ + (x - 1)];

}

}

}

if (y == new_height - 1) {

for (int x = 0; x

p[y * new_width + x] = p[(y - 2) * new_width + x];

}

for (int x = 0; x

p[x] = p[2 * new_width + x];

}

}

}

for (int y = 1; y

for (int x = 1; x

float value{ 0.f };

int count{ 0 };

for (int m = -1; m <= 1; ++m) {

for (int n = -1; n <= 1; ++n) {

value += p[(y + m) * new_width + (x + n)] * kernel[count++];

}

}

if (value

else if (value > 255.) dst_[(y - 1) * width_ + (x - 1)] = 255;

else dst_[(y - 1) * width_ + (x - 1)] = static_cast(value);

}

}

return 0;

}

} // namespace

int laplacian_cpu(const unsigned char* src, int width, int height, int ksize, unsigned char* dst, float* elapsed_time)

{

int ret{ -1 };

// ksize == 1: kernel={ 0, 1, 0, 1, -4, 1, 0, 1, 0 }

// ksize == 3: kernel={ 2, 0, 2, 0, -8, 0, 2, 0, 2 }

CHECK(ksize == 1 || ksize == 3);

//TIME_START_CPU

ret = Laplacian(src, dst, width, height, ksize);

//TIME_END_CPU

return ret;

}

main.cpp:#include "funset.hpp"

#include

#include

#include

#include

#include

#include

#include "common.hpp"

int test_image_process_laplacian()

{

cv::Mat src = cv::imread("E:/GitCode/CUDA_Test/test_data/images/lena.png", 0);

if (!src.data || src.channels() != 1) {

fprintf(stderr, "read image fail\n");

return -1;

}

int width{ 400 }, height{ 400 };

cv::resize(src, src, cv::Size(width, height));

std::unique_ptr data1(new unsigned char[width * height]), data2(new unsigned char[width * height]);

float elapsed_time1{ 0.f }, elapsed_time2{ 0.f }; // milliseconds

int ksize{ 1 };

CHECK(laplacian_cpu(src.data, width, height, ksize, data1.get(), &elapsed_time1) == 0);

//CHECK(laplacian_gpu(src.data, width, height, data2.get(), &elapsed_time2) == 0);

//fprintf(stdout, "gray image edge detection: laplacian: cpu run time: %f ms, gpu run time: %f ms\n", elapsed_time1, elapsed_time2);

cv::Mat dst;

cv::Laplacian(src, dst, src.depth(), ksize);

cv::imwrite("E:/GitCode/CUDA_Test/test_data/images/laplacian.png", dst);

CHECK(compare_result(data1.get(), dst.data, width*height) == 0);

//CHECK(compare_result(data1.get(), data2.get(), width*height) == 0);

save_image(src, dst, width, height / 2, "E:/GitCode/CUDA_Test/test_data/images/laplacian_result.png");

return 0;

}

执行结果如下：

由结果可知：C++实现结果和调用OpenCV的接口结果是完全一致的。