在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(精确到可忽略不计的数值差异)?

我正在研究一个特定的应用程序，在那里我需要实际的内核值，简单地找到另一种通过近似高斯差分来实现日志过滤的方法并不是我所追求的。

谢谢!

Answer 1

您可以使用公式

手动生成它

日志(x, y) =(1/(π*σ^ 4))* (1 - (x ^ 2 + y ^ 2)/(σ^ 2))* (e ^ (- (x ^ 2 + y ^ 2)/2σ^ 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

我要感谢老ufo把我推向正确的方向。我希望我不必通过做一个快速的matlab->openCV转换来重新发明轮子，但我猜这是我为快速解决方案提供的最佳解决方案。

注意-我只对方形内核这样做(否则很容易修改，但我不需要这样做…)。也许这篇文章可以用一种更优雅的形式写出来，但我做得很快，这样我就可以继续处理更紧迫的事情了。

From main function:

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过滤器内核来检测以像素为单位的预定义半径的blob。

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}{pisigma^4}left(1 - frac{x^2 + y^2}{2sigma^2}right)e^{frac{-(x^2+y^2)}{2sigma^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}{2sigma^2} = 0$
    # i.e. r^2 = 2sigma^2$ and thus $sigma = r / sqrt{2}$
    sigma = r_in_px/np.sqrt(2)
    # ksize such that filter covers $3sigma$
    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