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