2D Multi-scale dot enhancement filter based on Gaussian filter and Hessian matrix

Question

I'm trying to implement an algorithm that enhance dot and line like structures based on Hessian matrix, the algorithm uses Gaussian filter with different scales before calculating the hessian matrix for every pixel and later the Eigen values.

I want to use the algorithm to enhance pulmonary nodules (dot like structures) to extract them for later classification.

The Gaussian filter is used to reduce noise and preserve objects with specific scales (diameters in my case).

Here is the algorithm that I'm trying to implement:

A = imread('sliceX.png');

scales = [3, 5, 7, 11, 17]; % the scales in number of pixels 
% correspond to diameter of nodules in mm [1, 1.6, 2.4, 3.8, 6]
Ecircle = zeros(512, 512, length(scales)); % Ecircle will store the 
% results of each enhancement scale.

for sn=1:length(scales) % sn for scale number
% Smooth the original 2D image with a 2D Gaussian function of scale 
% Sigma_s.
    B = imgaussfilt(A, scales(sn));

    [gx, gy] = gradient(double(B));
    [gxx, gxy] = gradient(gx);
    [gxy, gyy] = gradient(gy); % it's normaly [gyx, gyy] but since 
    % gyx = gxy it's ok

    % loop through the image to calculte the eigen values for every 
    % pixel and based on that we choose the value of each scale 
    %enhancement filter
    for x = 1:512
        for y = 1:512
            % construct 2*2 hessian matrix for every pixel
            h = [gxx(x, y), gxy(x, y); gxy(x, y), gyy(x, y)]; 

            e = eig(h); % returns a vector with the two eigen values 
            % lambda1 and lambda2 in which lambda1 = e(1) and 
            % lambda2 = e(2)
            lambda1 = e(1); lambda2 = e(2);

            if abs(lambda1) < abs(lambda2)
                temp = lambda1;
                lambda1 = lambda2;
                lambda2 = temp;
            end

            if lambda1 < 0 && lambda2 < 0
                Ecircle(x, y, sn) = (abs(lambda2)^2)/abs(lambda1);
            end
        end
    end
    % multiply each enhancement scale by (sigma^2) as mentioned in the
    % article
    Ecircle(:, :, sn) = Ecircle(:, :, sn) * (scales(sn)^2);
end

I = max(Ecircle, [], 3);

This implementation for enhancing dot like structures (nodules), and can enhance line structures (vessels in my case) with just changing the if condition.

The problem is I'm not getting the supposed results that I should get:

The result of the developer of the method in the paper:

And this is my result: abvisioly the nodule in the image is greatly enhanced by there are great noise in the output image plus all the vessels joints are enhanced too (small and big)

I think the problem is with the Gaussian smoothing filter imgaussfilt(), and precisely with sigma or scales(sn) in this example. I've read that sigma should be in the same units of x and y i.e. number of pixels so I've transformed the diameters (that I got from the experimental results of the original article) to the number of pixels using PixelSpacing attribute form the original Dicom file metadata.

The diameters are [1, 1.6, 2.4, 3.8, 6] to cover approximately all the diameters of possible nodules.

Where I did get wrong and If the problem is with the value of Sigma,How can I specify sigma? and how can I use the imgaussfilt() function correctly?

Note: here is the image used in the example: input image

Answer 1

Questioner's answer...

1. sigma have the same units as x and y i.e. number of pixels.
2. In multi-scale filtering, the size of the filter must change when the sigma changes.
3. Obtain the number of pixels per one millimeter or the vice-versa. (I did this using the property of pixel spacing included in the DICOM metadata in Matlab you can do this as info=dicominfo('image.dcm'); and spacing=info.PixelSpacing;)

What I was doing wrong is not changing the filter size when I changed sigma in the different scales.

So, the solution to my problem is like this: the size of the gaussian kernel should be 2 times or preferably 3 times the value of sigma on eather side of the origin as Mr. Cris suggested. Gaussian filtering with Matlab's Image Processing Toolbox

cutoff = ceil(3*sigma);

In matlab there is two options, using fspecial() or imgaussfilt(), fspecial() is not recommended anymore in the newer versions of Matlab and latter is the recommended one.

h = fspecial('gaussian',2*cutoff+1,sigma);
B = conv2(A,h,'same');

or...

B = imgaussfilt(A, sigma, 'FilterSize', 2*cutoff+1);

Another thing that I did that corrected my results is converting my input image into type double using im2double() function.

Here is the result: