CODEWITHSUNDEEP
c++matlabopencvfilterkernel
CV_User
CV_User

Reputation: 1203

create 2D LoG kernel in openCV like fspecial in Matlab

My question is not how to filter an image using the laplacian of gaussian (basically using filter2D with the relevant kernel etc.).

What I want to know is how I generate the NxN kernel.

I'll give an example showing how I generated a [Winsize x WinSize] Gaussian kernel in openCV.

In Matlab:

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

In openCV:

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

Where sigma and WinSize are predefined.

I want to do the same for a Laplacian of Gaussian.

In Matlab:

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

How do I get the exact kernel in openCV (exact up to negligible numerical differences)?

I'm working on a specific application where I need the actual kernel values and simply finding another way of implementing LoG filtering by approximating Difference of gaussians is not what I'm after.

Thanks!

Upvotes: 4

Views: 7878

Answers (7)

check
check

Reputation: 13

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

Upvotes: 0

Saeed Masoomi
Saeed Masoomi

Reputation: 1834

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

Upvotes: 1

Erik Nijkamp
Erik Nijkamp

Reputation: 1

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

Upvotes: 0

user3667217
user3667217

Reputation: 2192

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

Upvotes: 0

br1
br1

Reputation: 443

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;

}

Upvotes: 1

CV_User
CV_User

Reputation: 1203

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

Upvotes: 1

old-ufo
old-ufo

Reputation: 2850

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

Upvotes: 4

Related Questions