added Hessian experimental

src/MATLAB/Hessian/drawEye.m

+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); % no mat2gray ! always quantize {x_i \in X} jointly 
+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;
+thickness = 8; % pixels
+sigma = 1/sqrt(3) * thickness * vx;
+
+% imLoG = HIP.LoG(im,sigma,0);
+plateRadii = 1/sqrt(3) * [0.5,8,8];
+% imPlate = HIP.LoG(im,plateRadii,0);
+% imGauss = HIP.Gaussian(im,sigma,1,0);
+% 
+if exist('v','var')
+    close(v)
+end
+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(:))]);
+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';   
+    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')
+    cb = colorbar();
+    cb.Limits = [cb1,cb2];
+%     cb.Ticks = [0,0.5,1];cb.TickLabels={'-1.0','0.0','1.0'};
+    f = getframe(gca);
+    v.writeVideo(f);
+    
+end
+close(v)
\ No newline at end of file

src/MATLAB/Hessian/goHessian.m

+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]);
+
+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 = imout .* imbinarize(imout);
+4;
\ No newline at end of file

src/MATLAB/Hessian/goHessian_2D.m

+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 = squeeze(im(600,:,:));
+im = im(:,:,100);
+im = mat2gray(im);
+voxelSize = CONSTANTS.imageData.PixelPhysicalSize(1:2);
+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,voxelSize);
+bw = imbinarize(im);
+bw = bwareaopen(bw,1e3);
+bw = imdilate(bw,ones(6));
+
+% im = mat2gray(im);
+Ig = HIP.Gaussian(im,0.5*[vx,0],1,0);
+[gx,gy] = gradient(Ig);
+[gxx,gxy] = gradient(gx);
+[gyx,gyy] = gradient(gy);
+
+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);gyx(idx1),gyy(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
+%         continue
+    end
+[x,y,z]=ind2sub(size(im),idx1);
+    ev = abs(ev);
+    rb =  ev(1) / ev(2) ;     
+    s = im(idx1);
+    v =  (exp(-rb^2 / 2)) * (1 - exp(-s^2 / 2));
+    imout(idx1) = v;
+end
+% imout = imout .* imbinarize(imout);
+4;
\ No newline at end of file

src/MATLAB/Hessian/loadImages.m

+function [im,imf,imn,imp,bw] = loadImages(fname,radii)
+
+im = leversc.loadImage(fname,1,1);
+% imf = RSF.imPreProcess(im,[3:8]);
+[imf,bw] = RSF.imPreProcess(im,radii);
+
+imp = imf;
+imp(imp<0) =0;
+imn=imf;
+imn(imf>0)=0;
+imn=abs(imn);
+
+