create 2D LoG kernel in openCV like fspecial in Matlab

Question 1

I want to thank old-ufo for nudging me in the correct direction. I was hoping I won't have to reinvent the wheel by doing a quick matlab-->openCV conversion but guess this is the best solution I have for a quick solution.

NOTE - I did this for square kernels only (easy to modify otherwise, but I have no need for that so...). Maybe this can be written in a more elegant form but is a quick job I did so I can carry on with more pressing matters.

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);

And the actual function is:

// 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;
}

Question 2

You can generate it manually, using formula

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));
     }

If function above is not OK, you can simply rewrite matlab version of 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

Question 3

There is some difference between your function and the matlab version: http://br1.einfach.org/tmp/log-matlab-vs-opencv.png.

Above is matlab fspecial('log', 31, 6) and below is the result of your function with the same parameters. Somehow the hat is more 'bent' - is this intended and what is the effect of this in later processing?

I can create a kernel very similar to the matlab one with these functions, which just directly reflect the LoG formula:

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;

}

Question 4

Here's a NumPy version that is directly translated from the fspecial function in MATLAB.

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()

Question 5

The code below is the exact equivalent to 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()

Question 6

I wrote exact Implementation of Matlab fspecial function in OpenCV

function:

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;
}

usage of this function :

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

compare OpenCV result with Matlab:

opencv result:

[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]

Matlab result for fspecial('gaussian', 4, 2.1) :

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

Question 7

Just for the sake of reference, here is a Python implementation which creates the LoG filter kernel to detect blobs of a pre-defined radius in pixels.

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