✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
In this paper we extend the Frangi filter[1] to recognize edges and do not enhance them. We give a theoretical framework for optimal scale selection and choice of the free parameters. We discuss discretization details concerning especially the discrete kernel used for building the scale-space and the choice of discrete scales. We present several experiments on phantom data to objectively and quantitatively compare and judge the filters. Experiments on real coronary angiograms enhance the improvement reached by the integration of the edge indicator.
⛄ 部分代码
function blobness = blobness3D(I, sigmas, spacing, tau, whiteondark)
% calculates blobness probability map (local sphericity) of a 3D input
% image
%
% blobness = blobness3D(V, sigmas, spacing, tau, whiteondark)
%
% inputs,
% I : 3D image
% sigmas : vector of scales on which the vesselness is computed
% spacing : input image spacing resolution - during hessian matrix
% computation, the gaussian filter kernel size in each dimension can
% be adjusted to account for different image spacing for different
% dimensions
% tau : (between 0.5 and 1) : parameter that controls response uniformity
% - lower tau -> more intense output response
% whiteondark: (true/false) : are vessels (tubular structures) bright on
% dark background or dark on bright (default for 3D is true)
%
% outputs,
% blobness: maximum blobness response over scales sigmas
%
% example:
% B = blobness3D(I, 1:5, [1;1;1], 1, true);
%
% Function was written by T. Jerman, University of Ljubljana (October 2014)
% Based on code by D. Kroon, University of Twente (May 2009)
verbose = 1;
if nargin<5
whiteondark = true; % default
end
I(~isfinite(I)) = 0;
I = single(I);
for j = 1:length(sigmas)
if verbose
disp(['Current Filter Sigma: ' num2str(sigmas(j)) ]);
end
[Lambda1, Lambda2, Lambda3] = volumeEigenvalues(I,sigmas(j),spacing,whiteondark);
% proposed filter
Lambda3M = Lambda3;
Lambda3M(Lambda3>=min(Lambda3(:)).*tau) = min(Lambda3(:)).*tau;
response = ((Lambda1.^2).*Lambda3M.* 27) ./ ((2*Lambda1+Lambda3M).^3);
response(abs(Lambda1)>abs(Lambda3M)) = 1;
response(Lambda1>=0) = 0;
response(Lambda2>=0) = 0;
response(Lambda3>=0) = 0;
response(~isfinite(response)) = 0;
%keep max response
if(j==1)
blobness = response;
else
blobness = max(blobness,response);
end
clear response Lambda1 Lambda2 Lambda3 Lambda3M
end
blobness = blobness ./ max(blobness(:));
blobness(blobness < 1e-2) = 0;
function [Lambda1, Lambda2, Lambda3] = volumeEigenvalues(V,sigma,spacing,whiteondark)
% calculates the three eigenvalues for each voxel in a volume
% Calculate 3D hessian
[Hxx, Hyy, Hzz, Hxy, Hxz, Hyz] = Hessian3D(V,sigma,spacing);
% Correct for scaling
c=sigma.^2;
Hxx = c*Hxx; Hxy = c*Hxy;
Hxz = c*Hxz; Hyy = c*Hyy;
Hyz = c*Hyz; Hzz = c*Hzz;
if whiteondark == false
c=-1;
Hxx = c*Hxx; Hxy = c*Hxy;
Hxz = c*Hxz; Hyy = c*Hyy;
Hyz = c*Hyz; Hzz = c*Hzz;
end
% reduce computation by computing vesselness only where needed
% S.-F. Yang and C.-H. Cheng, 揊ast computation of Hessian-based
% enhancement filters for medical images,?Comput. Meth. Prog. Bio., vol.
% 116, no. 3, pp. 215?25, 2014.
B1 = - (Hxx + Hyy + Hzz);
B2 = Hxx .* Hyy + Hxx .* Hzz + Hyy .* Hzz - Hxy .* Hxy - Hxz .* Hxz - Hyz .* Hyz;
B3 = Hxx .* Hyz .* Hyz + Hxy .* Hxy .* Hzz + Hxz .* Hyy .* Hxz - Hxx .* Hyy .* Hzz - Hxy .* Hyz .* Hxz - Hxz .* Hxy .* Hyz;
T = ones(size(B1));
T(B1<=0) = 0;
T(B2<=0) = 0;
T(B3<=0) = 0;
T(B1 .* B2 <= B3) = 0;
clear B1 B2 B3;
indeces = find(T==1);
Hxx = Hxx(indeces);
Hyy = Hyy(indeces);
Hzz = Hzz(indeces);
Hxz = Hxz(indeces);
Hyz = Hyz(indeces);
Hxy = Hxy(indeces);
% Calculate eigen values
[Lambda1i,Lambda2i,Lambda3i]=eig3volume(Hxx,Hxy,Hxz,Hyy,Hyz,Hzz);
% Free memory
clear Hxx Hyy Hzz Hxy Hxz Hyz;
Lambda1 = zeros(size(T));
Lambda2 = zeros(size(T));
Lambda3 = zeros(size(T));
Lambda1(indeces) = Lambda1i;
Lambda2(indeces) = Lambda2i;
Lambda3(indeces) = Lambda3i;
% some noise removal
Lambda1(~isfinite(Lambda1)) = 0;
Lambda2(~isfinite(Lambda2)) = 0;
Lambda3(~isfinite(Lambda3)) = 0;
Lambda1(abs(Lambda1) < 1e-4) = 0;
Lambda2(abs(Lambda2) < 1e-4) = 0;
Lambda3(abs(Lambda3) < 1e-4) = 0;
function [Dxx, Dyy, Dzz, Dxy, Dxz, Dyz] = Hessian3D(Volume,Sigma,spacing)
% This function Hessian3D filters the image with an Gaussian kernel
% followed by calculation of 2nd order gradients, which aprroximates the
% 2nd order derivatives of the image.
%
% [Dxx, Dyy, Dzz, Dxy, Dxz, Dyz] = Hessian3D(Volume,Sigma,spacing)
%
% inputs,
% I : The image volume, class preferable double or single
% Sigma : The sigma of the gaussian kernel used. If sigma is zero
% no gaussian filtering.
% spacing : input image spacing
%
% outputs,
% Dxx, Dyy, Dzz, Dxy, Dxz, Dyz: The 2nd derivatives
%
% Function is written by D.Kroon University of Twente (June 2009)
% defaults
if nargin < 2, Sigma = 1; end
if(Sigma>0)
%F=imbigaussian(Volume,Sigma,0.5);
F=imgaussian(Volume,Sigma,spacing);
else
F=Volume;
end
% Create first and second order diferentiations
Dz=gradient3(F,'z');
Dzz=(gradient3(Dz,'z'));
clear Dz;
Dy=gradient3(F,'y');
Dyy=(gradient3(Dy,'y'));
Dyz=(gradient3(Dy,'z'));
clear Dy;
Dx=gradient3(F,'x');
Dxx=(gradient3(Dx,'x'));
Dxy=(gradient3(Dx,'y'));
Dxz=(gradient3(Dx,'z'));
clear Dx;
function D = gradient3(F,option)
% This function does the same as the default matlab "gradient" function
% but with one direction at the time, less cpu and less memory usage.
%
% Example:
%
% Fx = gradient3(F,'x');
[k,l,m] = size(F);
D = zeros(size(F),class(F));
switch lower(option)
case 'x'
% Take forward differences on left and right edges
D(1,:,:) = (F(2,:,:) - F(1,:,:));
D(k,:,:) = (F(k,:,:) - F(k-1,:,:));
% Take centered differences on interior points
D(2:k-1,:,:) = (F(3:k,:,:)-F(1:k-2,:,:))/2;
case 'y'
D(:,1,:) = (F(:,2,:) - F(:,1,:));
D(:,l,:) = (F(:,l,:) - F(:,l-1,:));
D(:,2:l-1,:) = (F(:,3:l,:)-F(:,1:l-2,:))/2;
case 'z'
D(:,:,1) = (F(:,:,2) - F(:,:,1));
D(:,:,m) = (F(:,:,m) - F(:,:,m-1));
D(:,:,2:m-1) = (F(:,:,3:m)-F(:,:,1:m-2))/2;
otherwise
disp('Unknown option')
end
function I=imgaussian(I,sigma,spacing,siz)
% IMGAUSSIAN filters an 1D, 2D color/greyscale or 3D image with an
% Gaussian filter. This function uses for filtering IMFILTER or if
% compiled the fast mex code imgaussian.c . Instead of using a
% multidimensional gaussian kernel, it uses the fact that a Gaussian
% filter can be separated in 1D gaussian kernels.
%
% J=IMGAUSSIAN(I,SIGMA,SIZE)
%
% inputs,
% I: The 1D, 2D greyscale/color, or 3D input image with
% data type Single or Double
% SIGMA: The sigma used for the Gaussian kernel
% SIZE: Kernel size (single value) (default: sigma*6)
%
% outputs,
% J: The gaussian filtered image
%
% note, compile the code with: mex imgaussian.c -v
%
% example,
% I = im2double(imread('peppers.png'));
% figure, imshow(imgaussian(I,10));
%
% Function is written by D.Kroon University of Twente (September 2009)
if(~exist('siz','var')), siz=sigma*6; end
if(sigma>0)
% Filter each dimension with the 1D Gaussian kernels\
x=-ceil(siz/spacing(1)/2):ceil(siz/spacing(1)/2);
H = exp(-(x.^2/(2*(sigma/spacing(1))^2)));
H = H/sum(H(:));
Hx=reshape(H,[length(H) 1 1]);
x=-ceil(siz/spacing(2)/2):ceil(siz/spacing(2)/2);
H = exp(-(x.^2/(2*(sigma/spacing(2))^2)));
H = H/sum(H(:));
Hy=reshape(H,[1 length(H) 1]);
x=-ceil(siz/spacing(3)/2):ceil(siz/spacing(3)/2);
H = exp(-(x.^2/(2*(sigma/spacing(3))^2)));
H = H/sum(H(:));
Hz=reshape(H,[1 1 length(H)]);
I=imfilter(imfilter(imfilter(I,Hx, 'same' ,'replicate'),Hy, 'same' ,'replicate'),Hz, 'same' ,'replicate');
end
⛄ 运行结果
⛄ 参考文献
1. [T. Jerman, F. Pernus, B. Likar, Z. Spiclin, "*Enhancement of Vascular Structures in 3D and 2D Angiographic Images*", IEEE Transactions on Medical Imaging, 35(9), p. 2107-2118 (2016)
2. [T. Jerman, F. Pernus, B. Likar, Z. Spiclin, "*Blob Enhancement and Visualization for Improved Intracranial Aneurysm Detection*", IEEE Transactions on Visualization and Computer Graphics, 22(6), p. 1705-1717 (2016)
3. [T. Jerman, F. Pernus, B. Likar, Z. Spiclin, "*Beyond Frangi: an improved multiscale vesselness filter*", Proc. SPIE 9413, Medical Imaging 2015: Image Processing, 94132A (2015)
