在 openCV 中创建 2D LoG 内核，如 Matlab 中的 fspecial

Question

我的问题不是如何使用高斯的拉普拉斯算子过滤图像（基本上使用带有相关内核的 filter2D 等）。

我想知道的是我如何生成 NxN 内核。

我将举一个例子说明我是如何在 openCV 中生成 [Winsize x WinSize] 高斯内核的。

在 Matlab 中：

gaussianKernel = fspecial('gaussian', WinSize, sigma);

在openCV中：

cv::Mat gaussianKernel = cv::getGaussianKernel(WinSize, sigma, CV_64F);
cv::mulTransposed(gaussianKernel,gaussianKernel,false);

sigma 和 WinSize 是预定义的。

我想对高斯的拉普拉斯算子做同样的事情。

在 Matlab 中：

LoGKernel = fspecial('log', WinSize, sigma);

如何在 openCV 中获得精确的内核（精确到可以忽略的数值差异）？

我正在开发一个需要实际内核值的特定应用程序，并且只是通过近似高斯差来寻找另一种实现 LoG 过滤的方法，这不是我所追求的。

谢谢！

Answer 1

您可以使用公式手动生成它

LoG(x,y) = (1/(pi*sigma^4)) * (1 - (x^2+y^2)/(sigma^2))* (e ^ (- (x^ 2 + y^2) / 2sigma^2)

http://homepages.inf.ed.ac.uk/rbf/HIPR2/log.htm

cv::Mat kernel(WinSize,WinSize,CV_64F);
int rows = kernel.rows;
int cols = kernel.cols;
double halfSize = (double) WinSize / 2.0; 
for (size_t i=0; i<rows;i++)
  for (size_t j=0; j<cols;j++)
    { 
     double x = (double)j - halfSize;
     double y = (double)i - halfSize;
     kernel.at<double>(j,i) = (1.0 /(M_PI*pow(sigma,4))) * (1 - (x*x+y*y)/(sigma*sigma))* (pow(2.718281828, - (x*x + y*y) / 2*sigma*sigma));
     }

如果上面的功能不行，可以直接重写matlab版本的fspecial：

 case 'log' % Laplacian of Gaussian
 % first calculate Gaussian
 siz   = (p2-1)/2;
 std2   = p3^2;

 [x,y] = meshgrid(-siz(2):siz(2),-siz(1):siz(1));
 arg   = -(x.*x + y.*y)/(2*std2);

 h     = exp(arg);
 h(h<eps*max(h(:))) = 0;

 sumh = sum(h(:));
 if sumh ~= 0,
   h  = h/sumh;
 end;
 % now calculate Laplacian     
 h1 = h.*(x.*x + y.*y - 2*std2)/(std2^2);
 h     = h1 - sum(h1(:))/prod(p2); % make the filter sum to zero

Answer 2

我要感谢 old-ufo 将我推向正确的方向。 我希望我不必通过快速 matlab-->openCV 转换来重新发明轮子，但我想这是我快速解决方案的最佳解决方案。

注意 - 我只为方形内核执行此操作（否则很容易修改，但我不需要那样......）。 也许这可以写成更优雅的形式，但我做的很快，所以我可以继续处理更紧迫的事情。

来自主函数：

int WinSize(7); int sigma(1); // can be changed to other odd-sized WinSize and different sigma values
cv::Mat h = fspecialLoG(WinSize,sigma);

而实际作用是：

// return NxN (square kernel) of Laplacian of Gaussian as is returned by     Matlab's: fspecial(Winsize,sigma)
cv::Mat fspecialLoG(int WinSize, double sigma){
 // I wrote this only for square kernels as I have no need for kernels that aren't square
cv::Mat xx (WinSize,WinSize,CV_64F);
for (int i=0;i<WinSize;i++){
    for (int j=0;j<WinSize;j++){
        xx.at<double>(j,i) = (i-(WinSize-1)/2)*(i-(WinSize-1)/2);
    }
}
cv::Mat yy;
cv::transpose(xx,yy);
cv::Mat arg = -(xx+yy)/(2*pow(sigma,2));
cv::Mat h (WinSize,WinSize,CV_64F);
for (int i=0;i<WinSize;i++){
    for (int j=0;j<WinSize;j++){
        h.at<double>(j,i) = pow(exp(1),(arg.at<double>(j,i)));
    }
}
double minimalVal, maximalVal;
minMaxLoc(h, &minimalVal, &maximalVal);
cv::Mat tempMask = (h>DBL_EPSILON*maximalVal)/255;
tempMask.convertTo(tempMask,h.type());
cv::multiply(tempMask,h,h);

if (cv::sum(h)[0]!=0){h=h/cv::sum(h)[0];}

cv::Mat h1 = (xx+yy-2*(pow(sigma,2))/(pow(sigma,4));
cv::multiply(h,h1,h1);
h = h1 - cv::sum(h1)[0]/(WinSize*WinSize);
return h;
}

Answer 3

你的函数和matlab版本有一些区别： http://br1.einfach.org/tmp/log-matlab-vs-opencv.png。

上面是 matlab fspecial('log', 31, 6)，下面是具有相同参数的函数的结果。不知何故，这顶帽子更“弯曲”了——这是故意的吗？这在以后的处理中会产生什么影响？

我可以用这些函数创建一个与 matlab 非常相似的内核，它只是直接反映了 LoG 公式：

float LoG(int x, int y, float sigma) {
    float xy = (pow(x, 2) + pow(y, 2)) / (2 * pow(sigma, 2));
    return -1.0 / (M_PI * pow(sigma, 4)) * (1.0 - xy) * exp(-xy);
}

static Mat LOGkernel(int size, float sigma) {
   Mat kernel(size, size, CV_32F);
   int halfsize = size / 2;
   for (int x = -halfsize; x <= halfsize; ++x) {
        for (int y = -halfsize; y <= halfsize; ++y) {
            kernel.at<float>(x+halfsize,y+halfsize) = LoG(x, y, sigma);
        }
   }
   return kernel;

}

Answer 4

这是一个直接翻译自 MATLAB 中的fspecial 函数的 NumPy 版本。

import numpy as np
import sys


def get_log_kernel(siz, std):
    x = y = np.linspace(-siz, siz, 2*siz+1)
    x, y = np.meshgrid(x, y)
    arg = -(x**2 + y**2) / (2*std**2)
    h = np.exp(arg)
    h[h < sys.float_info.epsilon * h.max()] = 0
    h = h/h.sum() if h.sum() != 0 else h
    h1 = h*(x**2 + y**2 - 2*std**2) / (std**4)
    return h1 - h1.mean()

Answer 5

下面的代码完全等同于fspecial('log', p2, p3)：

def fspecial_log(p2, std):
   siz = int((p2-1)/2)
   x = y = np.linspace(-siz, siz, 2*siz+1)
   x, y = np.meshgrid(x, y)
   arg = -(x**2 + y**2) / (2*std**2)
   h = np.exp(arg)
   h[h < sys.float_info.epsilon * h.max()] = 0
   h = h/h.sum() if h.sum() != 0 else h
   h1 = h*(x**2 + y**2 - 2*std**2) / (std**4)
   return h1 - h1.mean()

Answer 6

我在 OpenCV 中编写了 Matlab fspecial 函数的精确实现

功能：

Mat C_fspecial_LOG(double* kernel_size,double sigma)
{
    double size[2]={  (kernel_size[0]-1)/2   , (kernel_size[1]-1)/2};
    double std = sigma;
    const double eps = 2.2204e-16;
    cv::Mat kernel(kernel_size[0],kernel_size[1],CV_64FC1,0.0);
    int row=0,col=0;
    for (double y = -size[0]; y <= size[0]; ++y,++row)
    {
       col=0;
       for (double x = -size[1]; x <= size[1]; ++x,++col)
       {
            kernel.at<double>(row,col)=exp( -( pow(x,2)  +  pow(y,2)  )  /(2*pow(std,2)));
       }
    }

    double MaxValue;
    cv::minMaxLoc(kernel,nullptr,&MaxValue,nullptr,nullptr);
    Mat condition=~(kernel < eps*MaxValue)/255;
    condition.convertTo(condition,CV_64FC1);
    kernel = kernel.mul(condition);

    cv::Scalar SUM = cv::sum(kernel);
    if(SUM[0]!=0)
    {
       kernel /= SUM[0];
    }

    return kernel;
} 

这个函数的用法：

double kernel_size[2] = {4,4};    // kernel size set to 4x4
double sigma =  2.1;
Mat kernel = C_fspecial_LOG(kernel_size,sigma);

将 OpenCV 结果与 Matlab 进行比较：

opencv 结果：

[0.04918466596701741, 0.06170341496034986, 0.06170341496034986, 0.04918466596701741;
 0.06170341496034986, 0.07740850411228289, 0.07740850411228289, 0.06170341496034986;
 0.06170341496034986, 0.07740850411228289, 0.07740850411228289, 0.06170341496034986;
 0.04918466596701741, 0.06170341496034986, 0.06170341496034986, 0.04918466596701741]

fspecial('gaussian', 4, 2.1) 的 Matlab 结果：

0.0492    0.0617    0.0617    0.0492
0.0617    0.0774    0.0774    0.0617
0.0617    0.0774    0.0774    0.0617
0.0492    0.0617    0.0617    0.0492

Answer 7

仅供参考，这里是一个 Python 实现，它创建 LoG 过滤器内核来检测以像素为单位的预定义半径的斑点。

def create_log_filter_kernel(r_in_px: float):
    """
    Creates a LoG filter-kernel to detect blobs of a given radius r_in_px.
    \[
        LoG(x,y) = \frac{-1}{\pi\sigma^4}\left(1 - \frac{x^2 + y^2}{2\sigma^2}\right)e^{\frac{-(x^2+y^2)}{2\sigma^2}}
    \]
    Look for maxima if blob is black, minima if blob is white.
    :param r_in_px:
    :return: filter kernel
    """
    # sigma from radius: LoG has zero-crossing at $1 - \frac{x^2 + y^2}{2\sigma^2} = 0$
    # i.e. r^2 = 2\sigma^2$ and thus $sigma = r / \sqrt{2}$
    sigma = r_in_px/np.sqrt(2)
    # ksize such that filter covers $3\sigma$
    ksize = int(np.round(sigma*3))*2 + 1
    # setup filter
    xgv = np.arange(0, ksize) - ksize / 2
    ygv = np.arange(0, ksize) - ksize / 2
    x, y = np.meshgrid(xgv, ygv)
    kernel = -1 / (np.pi * sigma**4) * (1 - (x**2 + y**2) / (2*sigma**2)) * np.exp(-(x**2 + y**2) / (2 * sigma**2))
    #normalize to sum zero (does not change zero crossing, I tried it out for r < 100)
    kernel -= np.sum(kernel) / ksize**2
    #this is important: normalize such that positive/negative parts are comparable over different scales
    kernel /= np.sum(kernel[kernel>0])
    return kernel