diff --git a/src/MATLAB/Hessian/drawEye.m b/src/MATLAB/Hessian/drawEye.m
new file mode 100644
index 0000000000000000000000000000000000000000..6621fcb8a6eb638ae19266b3b54b9b46d69faa57
--- /dev/null
+++ b/src/MATLAB/Hessian/drawEye.m
@@ -0,0 +1,46 @@
+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
diff --git a/src/MATLAB/Hessian/goHessian.m b/src/MATLAB/Hessian/goHessian.m
new file mode 100644
index 0000000000000000000000000000000000000000..c5b53a8f802ae5e55380c0efa6ddb70cc6c70c7f
--- /dev/null
+++ b/src/MATLAB/Hessian/goHessian.m
@@ -0,0 +1,71 @@
+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
diff --git a/src/MATLAB/Hessian/goHessian_2D.m b/src/MATLAB/Hessian/goHessian_2D.m
new file mode 100644
index 0000000000000000000000000000000000000000..b205641860ad02337ef78edba23a30d24c1ffc4c
--- /dev/null
+++ b/src/MATLAB/Hessian/goHessian_2D.m
@@ -0,0 +1,63 @@
+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
diff --git a/src/MATLAB/Hessian/loadImages.m b/src/MATLAB/Hessian/loadImages.m
new file mode 100644
index 0000000000000000000000000000000000000000..cc4480c4a2d925b9adcff94b64becab2e1a240b2
--- /dev/null
+++ b/src/MATLAB/Hessian/loadImages.m
@@ -0,0 +1,13 @@
+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);
+
+