clean up darktubes impreprocess

src/MATLAB/+RSF/darkTubes.m

@@ -3,7 +3,8 @@ function imfm = darkTubes(im,CONSTANTS,thickness)
 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;
-sigma = 1/sqrt(3) * thickness * vx;
+% radius pixels -> sigma LoG space
+sigma = 1/sqrt(3) * (thickness ./ 2) * vx;
 
 cudaTarget = getCudaTarget;
 imf = HIP.LoG(im,sigma,cudaTarget);
@@ -12,7 +13,7 @@ imn(imn>0) = 0;
 imn = abs(imn);
 imp=imf;imp(imp<0)=0;
 imfm = fibermetric(imp);
-imfm = imfm .* imregionalmax(imfm,6);
+imfm = imfm .* imregionalmax(imfm);
 % imfm = RSF.clipRSF(imfm);
 
 4;

src/MATLAB/+RSF/imPreProcess.m

 function [imf,bw] = imPreProcess(im,radii)
 
 % threshold mask
-bw = imbinarize(im);
-bw = bwareaopen(bw,1e3);
 if isempty(radii)
     % radii = {} ? return thresholded input image
+    bw = imbinarize(im);
+    bw = bwareaopen(bw,1e3);
     im(~bw) = 0;
     imf = im;
     return
@@ -17,16 +17,8 @@ cudaTarget = getCudaTarget;
 imf = HIP.LoG(im,sx,cudaTarget);
 [~,imp,imn] = LoG.separateLoG(imf);
 
-% threshold mask
-bw = imbinarize(im);
-bw = bwareaopen(bw,1e3);
-bw = imdilate(bw,ones(ceil(AX .* max(radii))));
-% ??? threshold before or after regionalmax ???
-% challenge problem -- connect this choice to K(X|M) 
 imn = imn .* imregionalmax(imn);
-imn = imn .* bw;
 imp = imp .* imregionalmax(imp);
-imp = imp .* bw;    
 
 imf = imp;
 imf(imn>=imp) = -1 * imn(imn>=imp);

src/MATLAB/+RSF/writeLever.m

+% e.g. %   
+leverFolder = '/g/leverjs/Schuman_OCT/viz'
+%   RSF.writeLever(imf{1},fname)
+% function writeLever(imf,fname)
+% load('/home/ac/git/rsf/src/MATLAB/clusterExamples/P10125_LR_8_qualYes.mat')
+% [imraw, CONSTANTS] = leversc.loadImage(fullfile(ROOT,tblMeta.filename{1}));
+% imf = RSF.darkTubes(imraw,CONSTANTS,[5]);
+
+
+
+% imImport = imraw;
+% 
+% imImport(:,:,:,2) = imf;
+% 
+% imImport = permute(imImport,[2,3,1,4]);
+imx = [];
+imx(:,:,:,1)=im;
+imx(:,:,:,2)=imn;
+imx(:,:,:,3)=imp;
+imx(:,:,:,4)=imfmX;
+imd = Import.importImage(imx,[]);
+
+CONSTANTS.imageData.DatasetName = ['imfm_' CONSTANTS.imageData.DatasetName];
+CONSTANTS.imageData.NumberOfChannels=4
+CONSTANTS.imageData.ChannelColors=[1,1,1;1,1,0;1,0,0;0,1,1];
+CONSTANTS.imageData.ChannelNames = {'raw';'LoG-';'LoG+';'fibermetric(LoG+) X'}
+
+[~,f,~] = fileparts(tblMeta.filename{1});
+Import.importImage(imx,CONSTANTS.imageData,leverFolder)
\ No newline at end of file

src/MATLAB/clusterExamples/getClipFrame.m

+% slice & dice
+function imRGB = getClipFrame(imIn,clipLimits,cm)
+
+if ~exist('cm','var')
+    cm = parula();
+end
+
+if 1 == clipLimits.dim
+   im = max(imIn(clipLimits.range,:,:),[],1);
+elseif 2 == clipLimits.dim
+   im = max(imIn(:,clipLimits.range,:),[],2);
+else
+    im = max(imIn(:,:,clipLimits.range),[],3);
+end
+im = squeeze(im);
+% im = mat2gray(im);
+[im] = gray2ind(im,size(cm,1));
+imRGB = ind2rgb(im,cm);
+% orient
+imRGB = flip(permute(imRGB,[2,1,3]),2);
+4;

src/MATLAB/clusterExamples/goVideoDarkTubes.m

+
+ROOT = '/g/leverjs/Schuman_OCT/OCT/qualifiedYes';
+PID = 'P10010'
+thickness = 5
+flist = dir(fullfile(ROOT,['/' PID '*.LEVER']));
+mx = regexp({flist.name},'P(?<PID>\d+)_(?<scanType>.+)_(?<date>\d+-\d+-\d+)_(?<time>\d+-\d+-\d+)_(?<eye>\w\w)_(?<scanID>.+?)_.*.LEVER','names');
+tblMeta = struct2table(vertcat(mx{:}));
+tblMeta.PID = str2double(tblMeta.PID);
+tblMeta.folder = {flist.folder}';
+tblMeta.filename = {flist.name}';
+tblMeta = tblMeta(contains(tblMeta.scanType,'Optic'),:);
+
+iid = 1;
+fname =  fullfile(tblMeta.folder{iid},tblMeta.filename{iid})
+[im,CONSTANTS] = leversc.loadImage(fname,1,1);
+[imraw,~] = leversc.loadImage(fname,1,2);
+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;
+sigma = 1/sqrt(3) * thickness./2 * vx;
+
+cudaTarget = getCudaTarget;
+imf = HIP.LoG(im,sigma,cudaTarget);
+imn = imf;
+imn(imn>0) = 0;
+imn = abs(imn);
+imp=imf;imp(imp<0)=0;
+imfm = fibermetric((imp));
+% imfmn = fibermetric(imn);
+% fig = figure('units','normalized','outerposition',[0 0 1 1])
+outfile = 'metricEmbeddings.tif';
+% v = VideoWriter('metricEmbeddings.avi','Uncompressed AVI');
+% v.open();
+%
+imfmX = imfm.*imregionalmax(imfm);
+% imp = imp .* imregionalmax(imp,6);
+% imn = imn .* imregionalmax(imn,6);
+
+sliceThickness = 100;
+for Z = sliceThickness: 1   :size(im,1)-sliceThickness
+    
+    clipLimits=[];
+    clipLimits.dim = 1;
+    clipLimits.range = [Z-sliceThickness:Z+sliceThickness];
+    clipLimits.range(clipLimits.range>size(im,1))=[];
+    clipLimits.range(clipLimits.range<1)=[];
+    
+    im1 = getClipFrame(im,clipLimits,gray());
+    im2 = getClipFrame(imn,clipLimits);
+    im3 = getClipFrame(imp,clipLimits);
+    im4 = getClipFrame(imfm,clipLimits);
+    im5 = getClipFrame(imfmn,clipLimits);
+    im6 = getClipFrame(imraw,clipLimits,gray());
+    im7 = getClipFrame(imfmX,clipLimits);
+
+    imOut = [ im1,im2, im7 ; im6, im3, im4   ];
+    if sliceThickness == Z
+        imwrite(imOut,outfile);
+    else
+        imwrite(imOut,outfile,'writemode','append');
+    end
+%     imagesc(imOut);hold on;set(gca,'XTick',[],'YTick',[]);axis equal
+%     xlim([1,400])
+%     ylim([1,400])
+%     text(10,10,['Z = ' num2str(Z-sliceThickness) ':' num2str(Z+sliceThickness)],'FontSize',10,'Color','w')
+%     text(200-5,10,['denoised input'],'FontSize',10,'Color','w','HorizontalAlignment','right')
+%     text(400-5,10,['bright blob (5\mum)'],'FontSize',10,'Color','w','HorizontalAlignment','right')
+%     text(200-5,210,['dark blob (5\mum)'],'FontSize',10,'Color','w','HorizontalAlignment','right')
+%     text(400-5,210,['dark tube'],'FontSize',10,'Color','w','HorizontalAlignment','right')
+%   
+%     writeVideo(v,f);
+
+    
+
+end
+% imwrite(f.cdata,[],outfile,"tiff",)
+% v.close()
+%
+% clipLimits = [];
+% clipLimits.dim = 3;
+% clipLimits.range = [70:110];
+% imH = getClipFrame(im,clipLimits,gray());
+% imHn = getClipFrame(imn,clipLimits);
+% imHp = getClipFrame(imp,clipLimits);
+% imHfm = getClipFrame(imfm,clipLimits);
+%
+% figure;imagesc([imH;imHn;imHp;imHfm])
+%
+%
+% imfm = imfm .* imregionalmax(imfm);
+% % imfm(~bw) = 0;
+%
+% imfm = RSF.clipRSF(imfm);
+%
+%
+%
+%
+%
+% 4;