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
The problem is I'm not getting the supposed results that I should get:
I think the problem is with the Gaussian smoothing filter
imgaussfilt(), and precisely with
scales(sn) in this example.
I've read that sigma should be in the same units of
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