added Hessian

src/MATLAB/Hessian/drawEye.m

@@ -11,10 +11,31 @@ thickness = 8; % pixels
 sigma = 1/sqrt(3) * thickness * vx;
 
 % imLoG = HIP.LoG(im,sigma,0);
-plateRadii = 1/sqrt(3) * [0.5,8,8];
+% % imcomplement :: dark tubes in im -> bright tubes in imLoG; 
+% % complement makes bright tubes dark again. same below. 
+% imLoG = mat2gray((imLoG));
+% 
+% plateRadii = 1/sqrt(3) * [0.5,8,8];
 % imPlate = HIP.LoG(im,plateRadii,0);
+% imPlate = mat2gray(imcomplement(imPlate));
+% 
 % imGauss = HIP.Gaussian(im,sigma,1,0);
+% imGauss = mat2gray(imcomplement(imGauss));
+% 
+% thickness = 3; % pixels
+% sigma = 1/sqrt(3) * thickness * vx;
+% bw = imbinarize(im);
+% bw = bwareaopen(bw,1e3);
+% bw = imdilate(bw,ones(4*thickness));
 % 
+% imHess = getHessian(im,bw,sigma);
+% imHessLoG = getHessian(imcomplement(imLoG),bw,sigma);
+% imHessLoG = mat2gray(imHessLoG);
+% imLoG = imLoG .* bw;
+% imPlate = imPlate .* bw;
+% imGauss = imGauss .* bw;
+% imHessLoG = imHessLoG .* bw;
+
 if exist('v','var')
     close(v)
 end
@@ -22,24 +43,26 @@ v = VideoWriter('gaussEye.mp4','MPEG-4');
 
 v.Quality = 100;
 v.open();
-cb1 = min([min(im(:)),min(imGauss(:)),min(imPlate(:)),min(imLoG(:))]);
-cb2 = max([max(im(:)),max(imGauss(:)),max(imPlate(:)),max(imLoG(:))]);
+
+imCompositeFM = 0 * squeeze(imfm(1,:,:));
+
 for i = 1:size(im,1)
-%     imc1 = ind2rgb(gray2ind(squeeze(im(i,:,:))',255),gray(255));
-%     imc2 = ind2rgb(gray2ind(squeeze(imGauss(i,:,:))',255),parula(255));
-%     imc3 = ind2rgb(gray2ind(squeeze(imLoG(i,:,:))',255),parula(255));
-%     im2 = [imc1,imc2,imc3];
-    im2 = squeeze([im(i,:,:),imGauss(i,:,:),imLoG(i,:,:),imPlate(i,:,:)]);
-    im2 = im2';   
+    if all(0==vertcat(bw(i,:,:)),'all'),continue,end
+    imCompositeFM = max(imCompositeFM,squeeze(imfm(i,:,:))');
+    im2 = squeeze([im(i,:,:),imPlate(i,:,:),imLoG(i,:,:)])';
+    im2 = [im2,imCompositeFM];
+%     im22 = squeeze([imLoG(i,:,:),imHessLoG(i,:,:),imHess(i,:,:)])';
+%     im2 = [im21;im22];   
     clf;imagesc(im2);axis image;axis off
     text(5,5,['Y = ' num2str(i)],'color','w')
-    text(205,5,'Gauss Blob','color','w')
-    text(405,5,'LoG Blob','color','w')
-    text(605,5,'LoG Plate','color','w')
+    text(205,5,'LoG bright plate','color','w')
+    text(405,5,'LoG dark blob','color','w')    
+    text(605,5,'fibermetric (LoG dark blob) composite','color','w')
     cb = colorbar();
-    cb.Limits = [cb1,cb2];
-%     cb.Ticks = [0,0.5,1];cb.TickLabels={'-1.0','0.0','1.0'};
-    f = getframe(gca);
+    cb.Limits = [0,1];
+%     cb.Limits = [cb1,cb2];
+%     cb.Ticks = [0,1];cb.TickLabels={'0.0','1.0'};
+    f = getframe(gcf);
     v.writeVideo(f);
     
 end

src/MATLAB/Hessian/getHessian.m

+function imout = getHessian(im,bw,sigma)
+
+% im = mat2gray(im);
+Ig = HIP.Gaussian(im,sigma,1,0);
+[gx,gy,gz] = gradient(Ig);
+[gxx,gxy,gxz] = gradient(gx);
+[gyx,gyy,gyz] = gradient(gy);
+[gzx,gzy,gzz] = gradient(gz);
+
+% vx = voxelSize;
+% for d1 = 1:length(vx)
+%     for d2 = 1:length(vx)
+%         v_xyz(d1,d2) = vx(d1) * vx(d2);
+%     end
+% end
+
+imout   = 0 * im;
+
+idx = find(bw);
+for i = 1:length(idx)
+    idx1 = idx(i);    
+       hessian1 = [ gxx(idx1), gxy(idx1), gxz(idx1);gyx(idx1),gyy(idx1),gyz(idx1);...
+        gzx(idx1),gzy(idx1),gzz(idx1)];
+    % divide by voxelSize for anisotropy
+%        hessian1 = hessian1 ./ v_xyz;
+
+    [evec,eigval] = eig(hessian1);
+    
+    eigval = diag(eigval);
+    
+    if all(0==eigval),continue,end
+    if any(isnan(eigval)),continue,end
+    [~,idxev] = sort(abs(eigval));
+    ev = eigval(idxev);
+%     ev = ev ./ norm(ev);    
+    if any(ev(2:3) > 0)
+        % dark tubes here, see Frangi table 1, \lamdba_{2,3} both H+
+        imout(idx1) = eps;
+        continue
+    end
+    ev = abs(ev);
+    rb =  ev(1) / sqrt(ev(2)*ev(3)) ; 
+    ra = ev(2)/ev(3); 
+    s = im(idx1);
+    v =  (1-exp(-ra^2 / 0.5) ) * (exp(-rb^2 / 0.5)) * (1 - exp(-s^2 / 0.5));
+    imout(idx1) = v;
+end
+
+imout = mat2gray(imout);
\ No newline at end of file

src/MATLAB/Hessian/goHessian.asv

+ROOT = 'C:\Users\ANDYR\OneDrive - Drexel University\teaching\ctia\images\OCT';
+fname = 'P10010_Optic Disc Cube 200x200_3-15-2017_10-26-59_OS_sn0975_cube_raw.LEVER';
+fname = fullfile(ROOT,fname);
+
+[im,CONSTANTS] = leversc.loadImage(fname,1,1);
+im = mat2gray(im);
+voxelSize = CONSTANTS.imageData.PixelPhysicalSize;
+voxelSize([1,2]) = voxelSize([2,1]); % swap (x,y,z) to matlab column major (y,x,z)
+vx = max(voxelSize)./voxelSize;
+% [im,imf,imn,imp,bw] = loadImages(fname,[16,3,3]);
+bw = imbinarize(im);
+bw = bwareaopen(bw,1e3);
+bw = imdilate(bw,ones(6));
+
+% im = mat2gray(im);
+thickness = 4;
+Ig = HIP.Gaussian(im,thickness * vx,1,0);
+[gx,gy,gz] = gradient(Ig);
+[gxx,gxy,gxz] = gradient(gx);
+[gyx,gyy,gyz] = gradient(gy);
+[gzx,gzy,gzz] = gradient(gz);
+
+for d1 = 1:length(voxelSize)
+    for d2 = 1:length(voxelSize)
+        v_xyz(d1,d2) = voxelSize(d1) * voxelSize(d2);
+    end
+end
+
+imout   = 0 * im;
+
+idx = find(bw);
+for i = 1:length(idx)
+    idx1 = idx(i);    
+       hessian1 = [ gxx(idx1), gxy(idx1), gxz(idx1);gyx(idx1),gyy(idx1),gyz(idx1);...
+        gzx(idx1),gzy(idx1),gzz(idx1)];
+    % divide by voxelSize for anisotropy
+       hessian1 = hessian1 ./ v_xyz;
+
+    [evec,eigval] = eig(hessian1);
+    
+    eigval = diag(eigval);
+    
+    if all(0==eigval),continue,end
+    if any(isnan(eigval)),continue,end
+    [~,idxev] = sort(abs(eigval));
+    ev = eigval(idxev);
+%     ev = ev ./ norm(ev);    
+    if ev(end)>0
+        % dark vessel
+%         continue
+    end
+    if ev(end)<0
+        % bright vessel
+        imout(idx1) = eps;
+        continue
+    end
+
+    ev = abs(ev);
+    rb =  ev(1) / sqrt(ev(2)*ev(3)) ; 
+    ra = ev(2)/ev(3); 
+    s = im(idx1);
+    v =  (1-exp(-ra^2 / 0.5) ) * (exp(-rb^2 / 0.5)) * (1 - exp(-s^2 / 0.5));
+    imout(idx1) = v;
+end
+% imout = imout .* imbinarize(imout);
+4;
\ No newline at end of file

src/MATLAB/Hessian/goHessian.m

@@ -10,62 +10,8 @@ vx = max(voxelSize)./voxelSize;
 % [im,imf,imn,imp,bw] = loadImages(fname,[16,3,3]);
 
 thickness = 8; % pixels
-
-bw = imbinarize(im);
-bw = bwareaopen(bw,1e3);
-bw = imdilate(bw,ones(4*thickness));
-
-% im = mat2gray(im);
-
 sigma = 1/sqrt(3) * thickness * vx;
-Ig = HIP.Gaussian(im,sigma,1,0);
-[gx,gy,gz] = gradient(Ig);
-[gxx,gxy,gxz] = gradient(gx);
-[gyx,gyy,gyz] = gradient(gy);
-[gzx,gzy,gzz] = gradient(gz);
-
-vx = voxelSize;
-for d1 = 1:length(vx)
-    for d2 = 1:length(vx)
-        v_xyz(d1,d2) = vx(d1) * vx(d2);
-    end
-end
-
-imout   = 0 * im;
-
-idx = find(bw);
-for i = 1:length(idx)
-    idx1 = idx(i);    
-       hessian1 = [ gxx(idx1), gxy(idx1), gxz(idx1);gyx(idx1),gyy(idx1),gyz(idx1);...
-        gzx(idx1),gzy(idx1),gzz(idx1)];
-    % divide by voxelSize for anisotropy
-%        hessian1 = hessian1 ./ v_xyz;
-
-    [evec,eigval] = eig(hessian1);
-    
-    eigval = diag(eigval);
-    
-    if all(0==eigval),continue,end
-    if any(isnan(eigval)),continue,end
-    [~,idxev] = sort(abs(eigval));
-    ev = eigval(idxev);
-%     ev = ev ./ norm(ev);    
-    if ev(end)>0
-        % dark vessel
-%         continue
-    end
-    if ev(end)<0
-        % bright vessel
-        imout(idx1) = eps;
-        continue
-    end
 
-    ev = abs(ev);
-    rb =  ev(1) / sqrt(ev(2)*ev(3)) ; 
-    ra = ev(2)/ev(3); 
-    s = im(idx1);
-    v =  (1-exp(-ra^2 / 0.5) ) * (exp(-rb^2 / 0.5)) * (1 - exp(-s^2 / 0.5));
-    imout(idx1) = v;
-end
+imout = getHessian(im,thickness,sigma);
 % imout = imout .* imbinarize(imout);
 4;
\ No newline at end of file