diff --git a/src/MATLAB/Hessian/drawEye.m b/src/MATLAB/Hessian/drawEye.m
index 6621fcb8a6eb638ae19266b3b54b9b46d69faa57..412856fb21a17605902bd4562358a72d08f94444 100644
--- a/src/MATLAB/Hessian/drawEye.m
+++ b/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
diff --git a/src/MATLAB/Hessian/getHessian.m b/src/MATLAB/Hessian/getHessian.m
new file mode 100644
index 0000000000000000000000000000000000000000..aceb1292008523df6da237cfd018024d4e7d9ae4
--- /dev/null
+++ b/src/MATLAB/Hessian/getHessian.m
@@ -0,0 +1,49 @@
+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
diff --git a/src/MATLAB/Hessian/goHessian.asv b/src/MATLAB/Hessian/goHessian.asv
new file mode 100644
index 0000000000000000000000000000000000000000..17dedc7dae07fda6cfaa3d93603d29deb2f08548
--- /dev/null
+++ b/src/MATLAB/Hessian/goHessian.asv
@@ -0,0 +1,66 @@
+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
diff --git a/src/MATLAB/Hessian/goHessian.m b/src/MATLAB/Hessian/goHessian.m
index c5b53a8f802ae5e55380c0efa6ddb70cc6c70c7f..04a4069522484b24447c3ead7556c644953d9171 100644
--- a/src/MATLAB/Hessian/goHessian.m
+++ b/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