Merge pull request #1 from zhengdali1990/master

test hessian and anisotropic filter

test/BlobFilter3D.m

@@ -0,0 +1,129 @@
+function [Iout,whatScale,Voutx,Vouty,Voutz]=BlobFilter3D(I,options)
+%
+%
+% Written by D.Kroon University of Twente (May 2009)
+
+% Constants vesselness function
+% 
+% Constants vesselness function
+%
+% edited by Zhengda Li Aug 2017 for blob detection
+
+defaultoptions = struct('ScaleRange', [1.5 3], 'ScaleRatio', 0.5,...
+    'Alpha', 1,'Gamma',1, ...
+    'verbose',true,'BlackWhite',true);
+
+% Process inputs
+if(~exist('options','var'))
+    options=defaultoptions; 
+else
+    tags = fieldnames(defaultoptions);
+    for i=1:length(tags)
+         if(~isfield(options,tags{i})),  options.(tags{i})=defaultoptions.(tags{i}); end
+    end
+    if(length(tags)~=length(fieldnames(options)))
+        warning('FrangiFilter3D:unknownoption','unknown options found');
+    end
+end
+
+% Use single or double for calculations
+if(~isa(I,'double')), I=single(I); end
+
+sigmas=options.ScaleRange(1):options.ScaleRatio:options.ScaleRange(2);
+sigmas = sort(sigmas, 'ascend');
+
+% Frangi filter for all sigmas
+for i = 1:length(sigmas)
+    % Show progress
+    if(options.verbose)
+        disp(['Current Frangi Filter Sigma: ' num2str(sigmas(i)) ]);
+    end
+
+    % Calculate 3D hessian
+    [Dxx, Dyy, Dzz, Dxy, Dxz, Dyz] = Hessian3D(I,sigmas(i));
+
+    if(sigmas(i)>0)
+        % Correct for scaling
+        c=(sigmas(i)^2);
+        Dxx = c*Dxx; Dxy = c*Dxy;
+        Dxz = c*Dxz; Dyy = c*Dyy;
+        Dyz = c*Dyz; Dzz = c*Dzz;
+    end
+
+    % Calculate eigen values
+    if(nargout>2)
+        [Lambda1,Lambda2,Lambda3,Vx,Vy,Vz]=eig3volume(Dxx,Dxy,Dxz,Dyy,Dyz,Dzz);
+    else
+        [Lambda1,Lambda2,Lambda3]=eig3volume(Dxx,Dxy,Dxz,Dyy,Dyz,Dzz);
+    end
+
+    % Free memory
+    clear Dxx Dyy  Dzz Dxy  Dxz Dyz;
+    if(options.BlackWhite)
+        filterMat=(Lambda2 > 0)|(Lambda3 > 0)|(Lambda1 > 0);
+        Lambda1(filterMat)=0;
+        Lambda2(filterMat)=0;        
+        Lambda3(filterMat)=0;
+        Lambda1=abs(Lambda1);
+        Lambda2=abs(Lambda2);
+        Lambda3=abs(Lambda3);
+
+    else
+        filterMat=(Lambda1 < 0)|(Lambda2 < 0)|(Lambda3 < 0);
+
+        Lambda1(filterMat)=0;
+        Lambda2(filterMat)=0;        
+        Lambda3(filterMat)=0;
+    end
+
+    % Calculate absolute values of eigen values
+    Rb=sqrt(Lambda2.*Lambda1)./Lambda3;
+%    Rc=Lambda2.*Lambda1.*Lambda3;
+    % Second order structureness. S = sqrt(sum(L^2[i])) met i =< D
+    S = sqrt(Lambda1.^2+Lambda2.^2+Lambda3.^2);
+    B = 2*options.Alpha^2; 
+    G = 2*options.Gamma^2;
+ %   T3 =    exp(-(C/LambdaAbs3.^2));
+    % Free memory
+
+    %Compute Vesselness function
+%    T1 =    exp(-(Ra.^2./A));    
+    T2 = 1-   exp(-(Rb.^2./B));
+    T0  = (1-exp(-S.^2./G));
+
+%    keyboard
+    % Free memory
+    clear S B C G Rb
+    %Compute Vesselness function
+    Voxel_data = T0.*T2;
+
+    % Free memory
+%    clear T0 T1 T2 T3;
+
+    clear T0 T1 T2
+    Voxel_data(~isfinite(Voxel_data))=0;
+
+
+    % Add result of this scale to output
+    if(i==1)
+        Iout=Voxel_data;
+        if(nargout>1)
+            whatScale = ones(size(I),class(Iout));
+        end
+        if(nargout>2)
+            Voutx=Vx; Vouty=Vy; Voutz=Vz;
+        end
+    else
+        if(nargout>1)
+            whatScale(Voxel_data>Iout)=i;
+        end
+        if(nargout>2)
+            Voutx(Voxel_data>Iout)=Vx(Voxel_data>Iout);
+            Vouty(Voxel_data>Iout)=Vy(Voxel_data>Iout);
+            Voutz(Voxel_data>Iout)=Vz(Voxel_data>Iout);
+        end
+        % Keep maximum filter response
+        Iout=max(Iout,Voxel_data);
+    end
+end
+

test/CalcBaseField.m

@@ -0,0 +1,34 @@
+function field=CalcBaseField(delta)
+%calculate the cloud of a single normal element
+eleLen=ceil(delta)*6+1;
+[x1,y1,z1]=meshgrid(-(eleLen-1)/2:1:(eleLen-1)/2,-(eleLen-1)/2:1:(eleLen-1)/2,-(eleLen-1)/2:1:(eleLen-1)/2);
+field=mvnpdf([x1(:),y1(:),z1(:)], [0,0,0],[delta,delta,delta/40]);
+field=reshape(field,eleLen,eleLen,eleLen);
+
+% 
+% c=-16*log(0.1)*(delta-1)/(pi*pi);
+% field=zeros(eleLen,eleLen,eleLen);
+% X=0;
+% Y=0;
+% Z=1;
+% for i=1:eleLen
+%     for j =1:eleLen
+%         for k = 1:eleLen
+%             x=i-(eleLen+1)/2;
+%             y=j-(eleLen+1)/2;
+%             z=k-(eleLen+1)/2;
+%             l=sqrt(x*x+y*y+z*z);
+%             if l==0
+%                 field(i,j,k)=200/delta/delta/delta;
+%             else
+%                 theta=acos(abs((z*Z+x*X+y*Y)/(sqrt(X*X+Y*Y+Z*Z))/l));
+%                 kappa=1/(sin(theta)+0.001);
+% %                 field(i,j,k)=exp(-(l*l+c*kappa*kappa*4+z*z/3)/(delta*delta))*2;
+%                 field(i,j,k)=exp(-(l*l+c*kappa*kappa*4)/(delta*delta))*2;                     
+%             end            
+%             
+%         end
+%     end
+% end
+
+field=field./(sum(sum(sum(field))));

test/CalcElementField.m

@@ -0,0 +1,20 @@
+function field=CalcElementField(delta,position,base)
+%calculate the cloud of any element
+eleLen=ceil(delta)*6+1;
+theta=position(1);
+gamma=position(2);
+R=[  cos(theta), -sin(theta), 0;cos(gamma)*sin(theta), ...
+     cos(gamma)*cos(theta), -sin(gamma);...
+     sin(theta)*sin(gamma), cos(theta)*sin(gamma),  cos(gamma)];
+ R=R^-1;
+ r=(eleLen-1)/2;
+[xx,yy,zz]=meshgrid(-r:r,-r:r,-r:r);
+coord=[xx(:),yy(:),zz(:)]*R;
+
+coord=coord+r+1;
+coord(coord<1)=1;
+coord(coord>eleLen)=eleLen;
+coord=round(coord(:,2))+round(coord(:,1)-1)*(eleLen)+round(coord(:,3)-1)*(eleLen*eleLen);
+base=base(:);
+
+field=reshape(base(coord),eleLen,eleLen,eleLen);

test/CalcVotingField.m

@@ -0,0 +1,50 @@
+% 
+function field=CalcVotingField(input, X, Y, Z)
+delta=8;
+steps=80;
+thresh=0.005;
+%matlabpool('open');
+[LenX,LenY,LenZ]=size(input);
+field=zeros(LenX,LenY,LenZ);
+data=cell(steps,steps);
+
+base=CalcBaseField(delta);
+stepSize=2*pi/(steps-1);
+
+for k=0:stepSize/2:pi
+    for j=0:stepSize:2*pi  
+        data{round(k/stepSize*2+1),round(j/stepSize+1)}=...
+            CalcElementField(delta,[j,k],base);
+    end
+end
+% ndata=cell(LenZ,1);
+% for i=1:LenZ
+%     ndata{i}=zeros(LenX,LenY,LenZ);
+% end
+
+cJ=atan(Y./(X+eps))+pi.*single(X<0)+pi*2*single(X>0&Y<0);
+%cJ(cJ<0)=cJ(cJ<0)+2*pi;
+cJ=round(cJ./stepSize)+1;
+cK=acos(sqrt(X.^2+Y.^2)./sqrt(X.^2+Y.^2+Z.^2));
+cK=round(cK./stepSize*2)+1;
+len=(size(data{1,1})-1)./2;
+lenX=len(1);lenY=len(2);lenZ=len(3);
+
+for k=1:LenZ
+    for j=1:LenY
+        for i=1:LenX
+            if(abs(input(i,j,k))>thresh)
+                x=i;y=j;z=k;
+                lx=min(x-1,lenX); rx=min(lenX,LenX-x);
+                ly=min(y-1,lenY); ry=min(lenY,LenY-y);
+                lz=min(z-1,lenZ); rz=min(lenZ,LenZ-z);               
+                field(x-lx:x+rx,y-ly:y+ry,z-lz:z+rz)=field(x-lx:x+rx,y-ly:y+ry,z-lz:z+rz)+...
+                    data{cK(i,j,k),cJ(i,j,k)}(lenX+1-lx:lenX+1+rx,lenY+1-ly:lenY+1+ry,lenZ+1-lz:lenZ+1+rz)...
+                    .*input(i,j,k);
+            end
+        end        
+    end
+    disp(['Caculationg plane ', num2str(k)])
+end
+
+%matlabpool('close');

test/FrangiFilter2D.m

@@ -0,0 +1,121 @@
+function [outIm,whatScale,Direction] = FrangiFilter2D(I, options)
+% This function FRANGIFILTER2D uses the eigenvectors of the Hessian to
+% compute the likeliness of an image region to vessels, according
+% to the method described by Frangi:2001 (Chapter 2).
+%
+% [J,Scale,Direction] = FrangiFilter2D(I, Options)
+%
+% inputs,
+%   I : The input image (vessel image)
+%   Options : Struct with input options,
+%       .FrangiScaleRange : The range of sigmas used, default [1 8]
+%       .FrangiScaleRatio : Step size between sigmas, default 2
+%       .FrangiBetaOne : Frangi correction constant, default 0.5
+%       .FrangiBetaTwo : Frangi correction constant, default 15
+%       .BlackWhite : Detect black ridges (default) set to true, for
+%                       white ridges set to false.
+%       .verbose : Show debug information, default true
+%
+% outputs,
+%   J : The vessel enhanced image (pixel is the maximum found in all scales)
+%   Scale : Matrix with the scales on which the maximum intensity 
+%           of every pixel is found
+%   Direction : Matrix with directions (angles) of pixels (from minor eigenvector)   
+%
+% Example,
+%   I=double(imread ('vessel.png'));
+%   Ivessel=FrangiFilter2D(I);
+%   figure,
+%   subplot(1,2,1), imshow(I,[]);
+%   subplot(1,2,2), imshow(Ivessel,[0 0.25]);
+%
+% Written by Marc Schrijver, 2/11/2001
+% Re-Written by D.Kroon University of Twente (May 2009)
+
+defaultoptions = struct('FrangiScaleRange', [1 3], 'FrangiScaleRatio', 0.5, 'FrangiBetaOne', 0.5, 'FrangiBetaTwo', 15, 'verbose',true,'BlackWhite',0);
+
+% Process inputs
+if(~exist('options','var')) 
+    options=defaultoptions; 
+else
+    tags = fieldnames(defaultoptions);
+    for i=1:length(tags)
+         if(~isfield(options,tags{i})),  options.(tags{i})=defaultoptions.(tags{i}); end
+    end
+    if(length(tags)~=length(fieldnames(options)))
+        warning('FrangiFilter2D:unknownoption','unknown options found');
+    end
+end
+
+sigmas=options.FrangiScaleRange(1):options.FrangiScaleRatio:options.FrangiScaleRange(2);
+sigmas = sort(sigmas, 'ascend');
+
+beta  = 2*options.FrangiBetaOne^2;
+c     = 2*options.FrangiBetaTwo^2;
+
+% Make matrices to store all filterd images
+ALLfiltered=zeros([size(I) length(sigmas)]);
+ALLangles=zeros([size(I) length(sigmas)]);
+
+% Frangi filter for all sigmas
+for i = 1:length(sigmas)
+    % Show progress
+    if(options.verbose)
+        disp(['Current Frangi Filter Sigma: ' num2str(sigmas(i)) ]);
+    end
+
+    % Make 2D hessian
+    [Dxx,Dxy,Dyy] = Hessian2D(I,sigmas(i));
+
+    % Correct for scale
+    Dxx = (sigmas(i)^2)*Dxx;
+    Dxy = (sigmas(i)^2)*Dxy;
+    Dyy = (sigmas(i)^2)*Dyy;
+
+    % Calculate (abs sorted) eigenvalues and vectors
+    [Lambda2,Lambda1,Ix,Iy]=eig2image(Dxx,Dxy,Dyy);
+
+    % Compute the direction of the minor eigenvector
+    angles = atan2(Ix,Iy);
+
+    % Compute some similarity measures
+    Lambda1(Lambda1==0) = eps;
+    Rb = (Lambda2./Lambda1).^2;
+    S2 = Lambda1.^2 + Lambda2.^2;
+
+    % Compute the output image
+    Ifiltered = exp(-Rb/beta) .*(ones(size(I))-exp(-S2/c));
+
+    % see pp. 45
+    if(options.BlackWhite)
+        Ifiltered(Lambda1<0)=0;
+    else
+        Ifiltered(Lambda1>0)=0;
+    end
+    % store the results in 3D matrices
+    ALLfiltered(:,:,i) = Ifiltered;
+    ALLangles(:,:,i) = angles;
+end
+
+% Return for every pixel the value of the scale(sigma) with the maximum 
+% output pixel value
+if length(sigmas) > 1
+    [outIm,whatScale] = max(ALLfiltered,[],3);
+    outIm = reshape(outIm,size(I));
+
+    if(nargout>1)
+        whatScale = reshape(whatScale,size(I));
+    end
+    if(nargout>2)
+        Direction = reshape(ALLangles((1:numel(I))'+(whatScale(:)-1)*numel(I)),size(I));
+    end
+else
+    outIm = reshape(ALLfiltered,size(I));
+    if(nargout>1)
+            whatScale = ones(size(I));
+    end
+    if(nargout>2)
+        Direction = reshape(ALLangles,size(I));
+    end
+end
+%outIm=outIm./max(outIm(:));

test/Hessian2D.m

@@ -0,0 +1,33 @@
+function [Dxx,Dxy,Dyy] = Hessian2D(I,Sigma)
+%  This function Hessian2 Filters the image with 2nd derivatives of a 
+%  Gaussian with parameter Sigma.
+% 
+% [Dxx,Dxy,Dyy] = Hessian2(I,Sigma);
+% 
+% inputs,
+%   I : The image, class preferable double or single
+%   Sigma : The sigma of the gaussian kernel used
+%
+% outputs,
+%   Dxx, Dxy, Dyy: The 2nd derivatives
+%
+% example,
+%   I = im2double(imread('moon.tif'));
+%   [Dxx,Dxy,Dyy] = Hessian2(I,2);
+%   figure, imshow(Dxx,[]);
+%
+% Function is written by D.Kroon University of Twente (June 2009)
+
+if nargin < 2, Sigma = 1; end
+
+% Make kernel coordinates
+[X,Y]   = ndgrid(-round(3*Sigma):round(3*Sigma));
+
+% Build the gaussian 2nd derivatives filters
+DGaussxx = 1/(2*pi*Sigma^4) * (X.^2/Sigma^2 - 1) .* exp(-(X.^2 + Y.^2)/(2*Sigma^2));
+DGaussxy = 1/(2*pi*Sigma^6) * (X .* Y)           .* exp(-(X.^2 + Y.^2)/(2*Sigma^2));
+DGaussyy = DGaussxx';
+
+Dxx = imfilter(I,DGaussxx,'conv');
+Dxy = imfilter(I,DGaussxy,'conv');
+Dyy = imfilter(I,DGaussyy,'conv');