darktubes combined with imPreProcess

src/MATLAB/+RSF/LoG

+
+function [imf,imp,imn] = LoG(im,radii)
+
+cudaTarget = getCudaTarget;
+imf = HIP.LoG(im,sigma,cudaTarget);
+imn = imf;
+imn(imn>0) = 0;
+imn = abs(imn);
+imp=imf;imp(imp<0)=0;
+
+

src/MATLAB/+RSF/darkTubes.m

-function imOut = darkTubes(im,CONSTANTS,thicknessIn)
-
-thickness = abs(thicknessIn);
-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;
-% radius pixels -> sigma LoG space
-sigma = 1/sqrt(3) * (thickness ./ 2) * vx;
-
-cudaTarget = getCudaTarget;
-
-if thicknessIn < 0
-    % bright plate
-    sigmaPlate = [8,0.25,8]; %[1,10,1].*sigma;
-    imfPlate = HIP.LoG(im,sigmaPlate,cudaTarget);
-    imnPlate = imfPlate;
-    imnPlate(imnPlate>0) = 0;
-    imnPlate = abs(imnPlate);
-%     imnPlate = imnPlate .* imbinarize(imnPlate);
-%     imnPlate = mat2gray(imnPlate);
-    imnPlate = imnPlate .* imregionalmax(imnPlate);
-    imOut = imnPlate;
-else        
-    % dark tube
-    im = mat2gray(im);
-    imf = HIP.LoG(im,sigma,cudaTarget);
-    imp=imf;imp(imp<0)=0;
-    imfm = fibermetric(imp);
-    imfm = imfm .* imregionalmax(imfm);
-    imOut = imfm;
-end
-
-4;
-
-

src/MATLAB/clusterExamples/csf_left_right.m

+function csf = csf_left_right(tblMeta,A,idxPerm,idxGT)
+
+% keep idxPerm elements from tblMeta and corresponding distance matrix A
+if ~issymmetric(A) || ~all(diag(A))==0
+    warning('input A not symmetric');
+end
+
+tm = tblMeta(idxPerm,:);
+A = A(idxPerm,idxPerm);
+idxGT = idxGT(idxPerm);
+[~,Y] = Cluster.SpectralCluster(A,2);
+csf = CSF.csf_spatial(A,idxGT);
+
+% idxOD = find(strcmp(tm.eye,'OD'));
+% idxOS = find(strcmp(tm.eye,'OS'));
+% 
+% % draw scatter plot
+% figure(1)
+% clf;hold on
+% plot3(Y(idxOD,1),Y(idxOD,2),tm.dx(idxOD),'r*')
+% plot3(Y(idxOS,1),Y(idxOS,2),tm.dx(idxOS),'og')
+% lidx = {};
+% for i = 1:height(tm),lidx{i}=i;end
+% text(Y(:,1),Y(:,2),tm.dx(:),lidx)
+% zlabel('time')
+% legend({'OD','OS'})
+% xlabel('NCV1')
+% ylabel('NCV2')  
+% zlabel('time (months)')
+% title(['PID = ' num2str(tblMeta.PID(1)) ' :: normalized compression vectors idxPerm = '  jsonencode(idxPerm) ')'])
+% % compute CSF
+
+
+

src/MATLAB/clusterExamples/getClipFrame.m

@@ -5,17 +5,26 @@ if ~exist('cm','var')
     cm = parula();
 end
 
+clipLimits.range = max(clipLimits.range,ones(size(clipLimits.range)));
+im = mat2gray(imIn);
 if 1 == clipLimits.dim
-   im = max(imIn(clipLimits.range,:,:),[],1);
+   im = max(im(clipLimits.range,:,:),[],1);
 elseif 2 == clipLimits.dim
-   im = max(imIn(:,clipLimits.range,:),[],2);
+   im = max(im(:,clipLimits.range,:),[],2);
 else
-    im = max(imIn(:,:,clipLimits.range),[],3);
+    im = max(im(:,:,clipLimits.range),[],3);
 end
 im = squeeze(im);
-% im = mat2gray(im);
-[im] = gray2ind(im,size(cm,1));
-imRGB = ind2rgb(im,cm);
+STRIDE = min(size(imIn));
+% im clipped to 2-d
+if any([STRIDE,STRIDE] ~= size(im))
+    im = imresize(im,[STRIDE,STRIDE]); % careful here!
+end
+cmlength = length(cm);
+imInd = gray2ind(im,cmlength);
+imRGB = ind2rgb(imInd,cm);
 % orient
+% ACK verify
 imRGB = flip(permute(imRGB,[2,1,3]),2);
+
 4;

src/MATLAB/clusterExamples/get_pro.m

+function [imRGB,imInd] = get_pro(im,cmfn)
+       
+     
+    if ~exist('cmfn','var')
+        cmfn = @gray;
+    end
+    impro = squeeze(max(im,[],1));
+%     if max(impro(:))<0.1
+%         impro = mat2gray(impro);
+%     end
+%     
+impro = mat2gray(impro);
+    cmlength = min(length(unique(im)),256);
+    imInd = grayslice(impro,cmlength);
+    
+    imRGB = ind2rgb(imInd,cmfn(cmlength));
+    % orient
+    % ACK verify
+    imRGB = flip(permute(imRGB,[2,1,3]),2);
+    imInd = flip(permute(imInd,[2,1]),2);
+    imInd = imInd ./ 256;
+4;
+    
+
+

src/MATLAB/clusterExamples/goDrawDarkTubes.m

@@ -9,10 +9,12 @@ 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})
+
+% start here with fname path to .LEVER file
 [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)
@@ -25,47 +27,27 @@ imn = imf;
 imn(imn>0) = 0;
 imn = abs(imn);
 imp=imf;imp(imp<0)=0;
-% imp=imf;imp(imp<imn)=0;
-% imfb = mat2gray(imbinarize(imp));
 imfm = fibermetric((imp));
-% imfm = fibermetric(mat2gray(imp));
-% NOTE as with imPreProcess, regional max first, then threshold
-% imfm = imp;
-
-close all
-
-figure(1)
-
-clipLimits=[];
-clipLimits.dim = 1;
-clipLimits.range = [400:600];
-im1 = getClipFrame(im,clipLimits,gray());
-im2 = getClipFrame(imn,clipLimits);
-im3 = getClipFrame(imp,clipLimits);
-im4 = getClipFrame(imfm,clipLimits);
-
-imOut = [im1,im2;im3,im4];
-figure;axis equal;imagesc(imOut)
-set(gca,'XTick',[],'YTick',[])
-
-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);
+imfmX = imfm.*imregionalmax(imfm);
 
+% bright plate
+sigmaPlate = [5,1,1].*sigma;
+imfPlate = HIP.LoG(imn,sigmaPlate,cudaTarget);
+imnPlate = imfPlate;
+imnPlate(imnPlate>0) = 0;
+imnPlate = abs(imnPlate);
+imnPlateX = imnPlate.*imregionalmax(imnPlate,6);
 
+im1 = get_pro(im);
+imn1 = get_pro((imn));
+imp1 = get_pro((imp));
+imfm1 = get_pro((imfm),@parula);
+[imfmX1,imfmX1_ind] = get_pro(imfmX,@parula);
+imraw1 = get_pro((imraw));
+[implate1,implate1_ind] = get_pro(imnPlate,@parula);
+[imtube1,imfm1_ind] = get_pro(imfm,@parula);
+[~,imXplate1] = get_pro(imnPlateX,@parula);
 
+imOut = [imraw1,imn1, implate1; imp1,imtube1,imfmX1];
 
 
-4;

src/MATLAB/clusterExamples/goLR.m

 tic
 ROOT = '/g/leverjs/Schuman_OCT/OCT/qualifiedYes';
 % ROOT = '/g/leverjs/Schuman_OCT/OCT/qualified';
-
-PID = 'P10010';
+% ROOT = '/g/leverjs/Schuman_OCT/OCT/loanCX/qualified';
+% ROOT = '/g/leverjs/Schuman_OCT/OCT/combined';
+% PID = 'P10010';
 % PID = 'P10125'
+PID = 'P10091'
+% PID = 'P10053'
 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{:}));
@@ -21,8 +24,9 @@ p = ljsStartParallel(96);
 % radii = {[0.5,8,8]}; % plate
 % radii = {[5,5,5]};
 
-radii = {[10]} % dark tubes
-% radii = {[0.5,8,8],[5,5,5],[8]}; % plate,blob, dark tube
+% radii = {[12]} % dark tubes
+% radii = {[10],[0.5,8,8]}; % plate,blob
+radii = {[5]}; % plate,blob, dark tube
 % radii = {[5,5,5],[8]}; % plate,blob, dark tube
 radii
 
@@ -40,6 +44,7 @@ end
 
 A = Cluster.Regularize(d);
 [~,Y] = Cluster.SpectralCluster(A,2);
+[Y,e] = cmdscale(A)
 toc
 
 % compute time for scatter plot 3rd axis
@@ -54,6 +59,9 @@ figure(1)
 clf;hold on
 plot3(Y(idxOD,1),Y(idxOD,2),tblMeta.dx(idxOD),'r*')
 plot3(Y(idxOS,1),Y(idxOS,2),tblMeta.dx(idxOS),'og')
+lidx = {};
+for i = 1:height(tblMeta),lidx{i}=i;end
+text(Y(:,1),Y(:,2),tblMeta.dx(:),lidx)
 zlabel('time')
 legend({'OD','OS'})
 xlabel('NCV1')

src/MATLAB/clusterExamples/goSelectLR.m

+
+csf2 = [];
+qq = readtable('/home/ac/git/rsf/src/MATLAB/qualify/nyu onh matched.xlsx');
+qq = qq(find(strcmp(qq.PROJECT_ID,'10091')),:);
+
+for i = 1:height(tblMeta)
+    idx = find(cellfun(@length,strfind(qq.file_name,tblMeta.filename{i}(1:end-6))));
+    if isempty(idx)
+        tblMeta.qs{i} = '';
+    else
+        tblMeta.qs{i} = qq.qualification_status{idx};
+    end
+end
+    
+for idxDel = 1:size(A,1)
+    idxPerm = [1:size(A,1)];
+    idxPerm(idxDel) = [];
+    csf = csf_left_right(tblMeta,A,idxPerm,idxGT);
+    csf2(idxDel) = csf(2)
+    4;
+end
+
+A = Cluster.Regularize(d);
+[~,Y] = Cluster.SpectralCluster(A,2);
+toc
+
+% compute time for scatter plot 3rd axis
+tblMeta.dx =    datetime(tblMeta.date,'InputFormat','MM-dd-yyyy');
+tblMeta.dx = calmonths(between(min(tblMeta.dx),tblMeta.dx));
+
+idxOD = find(strcmp(tblMeta.eye,'OD'));
+idxOS = find(strcmp(tblMeta.eye,'OS'));
+
+% draw scatter plot
+figure(1)
+clf;hold on
+plot3(Y(idxOD,1),Y(idxOD,2),tblMeta.dx(idxOD),'r*')
+plot3(Y(idxOS,1),Y(idxOS,2),tblMeta.dx(idxOS),'og')
+lidx = {};
+for i = 1:height(tblMeta),lidx{i}=[num2str(i) '_' tblMeta.qs{i}];end
+text(Y(:,1),Y(:,2),tblMeta.dx(:),lidx,'Interpreter','none','FontSize',8)
+zlabel('time')
+legend({'OD','OS'})
+xlabel('NCV1')
+ylabel('NCV2')  
+zlabel('time (months)')
+% title(['PID = ' PID ' :: normalized compression vectors (radii = '  jsonencode(radii) ')'])
\ No newline at end of file

src/MATLAB/ncdVSvfmd/embedSameDate.m

+dNCD = [];
+
+mxx = dxx{id_dxx};
+idMin = min(mxx(:,2));
+
+mxx(:,2) = mxx(:,2) - idMin + 1;
+mxx(:,3) = mxx(:,3) - idMin + 1;
+d=[];
+for i = 1:size(dxx{id_dxx},1)
+        d(mxx(i,2),mxx(i,3)) = mxx(i,1);
+end
+
+A = Cluster.Regularize(d);
+
+[idx,Y] = Cluster.SpectralCluster(A,size(A,1));
+
+tx = pData.dx(unique(mxx(:,2))+idMin-1);
+
+tx = calmonths(between(min(tx),tx));
+
+mdl = fitlm(Y,tx,'linear');
+txPred = predict(mdl,Y);
+ePrediction = abs(tx-txPred);
+
+

src/MATLAB/ncdVSvfmd/getRNFL.m

 function iRNFL = getRNFL(p1,tblRNFL)
 
+return -1;
+
 % find p1 in tblRNFL and return avg thickness
 [~,f1,~] = fileparts(p1.filename);
 [~,tblFname,~] = fileparts(tblRNFL.file_name);

src/MATLAB/ncdVSvfmd/getStats.m

 
-function tblVelo = getStats(dxx, pData, tblMD, tblRNFL)
-
-if ~exist('tblRNFL','var')
-    tblRNFL = [];
-end
+function tblVelo = getStats(dxx, pData, tblMD)
 
 tblVelo = table();
 idxUnused = 1:height(pData);
 for i = 1:length(dxx)
     
     for j = 1:size(dxx{i},1)
-        jRNFL = getRNFL(pData(j,:),tblRNFL);
         m = dxx{i}(j,2);
         n = dxx{i}(j,3);
         % m and n are indices into pData
@@ -33,10 +28,6 @@ for i = 1:length(dxx)
             continue
         end
 
-        mRNFL = getRNFL(pData(m,:),tblRNFL);
-        nRNFL = getRNFL(pData(n,:),tblRNFL);
-        deltaRNFL = nRNFL - mRNFL;
-
         nt = table();
         nt.subjectID = pData.subjectID(m);
         idxm = find(dxx{i}(:,2)==m & dxx{i}(:,3) == m);
@@ -48,9 +39,6 @@ for i = 1:length(dxx)
         nt.mn = [m,n];                
         nt.ddx = pData.dx(n) - pData.dx(m);
         
-        nt.mRNFL = mRNFL;
-        nt.nRNFL = nRNFL;
-        nt.deltaRNFL = deltaRNFL;
         % find sex
         idx = find(tblMD.subjectID == pData.subjectID(m),1,'first');
         nt.sex = tblMD.gender_full(idx);

src/MATLAB/ncdVSvfmd/goKernelCorr.m

@@ -16,8 +16,8 @@ tblMeta = getMetaOCT(ROOT,eye,scanType);
 % load('../qualify/onh_rnfl_thickness.mat');
 
 % radii = {[77,5,5]}; % 
-radii = {[5]}; % 
-% radii = {[0.5,8,8]}
+% radii = {[5]}; % 
+radii = {-1*[8,0.5,8]}
 radii
 
 [rho2,tblVelo, pData, dxx]  = kernelCorr(radii,tblVFMD,tblMeta);
@@ -25,5 +25,5 @@ radii
 % optional train model here, e.g.
 % [rmse1,~] = regnetCorr(tblVelo);median(rmse1)
 
-% saveFile = fullfile(SAVE_FOLDER,['qualEns_ncdVSvfmd_' jsonencode(radii) '.mat']);
+% saveFile = fullfile(SAVE_FOLDER,['qualYes_ncdVSvfmd_' jsonencode(radii) '.mat']);
 % save(saveFile)

src/MATLAB/ncdVSvfmd/kernelCorr.m

-function [rho2,tblVelo, pData, dxx] = kernelCorr(radii,tblVFMD,tblMeta,tblRNFL)
+function [rho2,tblVelo, pData, dxx] = kernelCorr(radii,tblVFMD,tblMeta)
 
 kernelCorrStart = tic();
 
@@ -30,7 +30,8 @@ pData = table();
 for i = 1:length(mdInfo)
     for pj = 1:length(mdInfo(i).dx)
         idx = find(tblMeta.subjectID == mdInfo(i).id & tblMeta.dx == mdInfo(i).dx(pj));
-        p1 = tblMeta(idx,:);
+%         p1 = tblMeta(idx(1),:); % NOTE - AC choose at random 1 from same date stacks
+        p1 = tblMeta(idx,:); % NOTE - AC choose at random 1 from same date stacks
         p1.md = repmat(mdInfo(i).md(pj),height(p1),1);
         pData = [pData;p1];
     end
@@ -85,8 +86,8 @@ parfor wx = 1:length(workList)
     end
 end
 
-tblVelo = getStats(dxx,pData,tblVFMD,tblRNFL);
-[rho,pCorr] = corr((tblVelo.ncd),(tblVelo.dmd));
+tblVelo = getStats(dxx,pData,tblVFMD);
+[rho,pCorr] = corr((tblVelo.ncd),(tblVelo.dmd),'type','Spearman');
 rho2 = rho^2
 toc(kernelCorrStart)
 4;

src/MATLAB/nonHumanPrimate/goHuman.m

-tic
-path(path,fullfile('../clusterExamples/'))
-% ROOT = '/g/leverjs/Schuman_OCT/OCT/nonHumanPrimate';
-% fname = fullfile(ROOT,'31344_OS_V_5x5_0_0004257.LEVER');
-% ROOT = '/g/leverjs/Schuman_OCT/OCT/nonHumanPrimate';
-% fname = fullfile(ROOT,'31344_OS_V_5x5_0_0004257.LEVER');
-
-
-
-[imraw,CONSTANTS] = leversc.loadImage(fname,1,2);
-% imraw = flip(permute(imraw,[2,1]),2); % ACK verify direction!
-rx = 5;
-im_medfilt = medfilt3(imraw,[5,5,5]);
-% note NLM takes minutes (vs. e.g. seconds for LoG etc)...use batch!
-cudaTarget = getCudaTarget;
-im = HIP.NLMeans(im_medfilt,0.1,12,1,cudaTarget);
-
-thickness = 5;
-
-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);
-imp=imf;imp(imp<0)=0;
-
-imn = imf;
-imn(imn>0)=0;imn=abs(imn);
-
-sigmaPlate = [1,5,1].*sigma;
-imfPlate = HIP.LoG(imn,sigmaPlate,cudaTarget);
-imnPlate = imfPlate;
-imnPlate(imnPlate>0) = 0;
-imnPlate = abs(imnPlate);
-imnPlate = mat2gray(imnPlate);
-
-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));
-
-outfile = 'nhp_metricEmbeddings.tif';
-
-%
-imfmX = imfm.*imregionalmax(imfm);
-
-
-
-% bright plate
-sigmaPlate = [1,5,1].*sigma;
-imfPlate = HIP.LoG(im,sigmaPlate,cudaTarget);
-
-imnPlate = imfPlate;
-imnPlate(imnPlate>0) = 0;
-imnPlate = abs(imnPlate);
-imnPlateX = imnPlate.*imregionalmax(imnPlate,6);
-imnPlateX = mat2gray(imnPlateX);
-% impPlate = imfPlate;
-% impPlate(imnPlate<0) = 0;
-% impPlateX = impPlate.*imregionalmax(imnPlate,6);
-
-sliceThickness = 10;
-bClear = true;
-for Z = 100: 1   :size(im,1)-100
-    
-    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);    
-    im6 = getClipFrame(imraw,clipLimits,gray());
-%     im7 = getClipFrame(imfmX,clipLimits);
-    
-    im7 = getClipFrame(imnPlateX,clipLimits);
-    im7 = imdilate(im7,strel('disk',3));
-    im8 = getClipFrame(imnPlate,clipLimits);
-
-    imOut = [ im6,im2, im8 ; im3, im4, im7   ];
-    clf;imagesc(imOut);
-    title(num2str(Z))
-    drawnow
-    if bClear
-        imwrite(imOut,outfile);
-        bClear = false;
-    else
-        imwrite(imOut,outfile,'writemode','append');
-    end
-end
-toc
\ No newline at end of file

src/MATLAB/nonHumanPrimate/renderPipeline.m

@@ -15,9 +15,9 @@ for Z = 1 : 1   : size(im,1)
     imPipeline = [...
         getClipFrame(imraw,clipLimits,gray()),...
         getClipFrame(imn,clipLimits),...
-        getClipFrame(imp,clipLimits),...
         getClipFrame(imnPlate,clipLimits),...  
-        getClipFrame(imfmX,clipLimits),...
+        getClipFrame(imp,clipLimits),...
+        getClipFrame(imfm,clipLimits),...
         ];    
     
     for idim = 2:3        
@@ -27,9 +27,9 @@ for Z = 1 : 1   : size(im,1)
         imPipeline = [imPipeline;...
             getClipFrame(imraw,clipLimits,gray()),...
             getClipFrame(imn,clipLimits),...
-            getClipFrame(imp,clipLimits),...
             getClipFrame(imnPlate,clipLimits),...
-            getClipFrame(imfmX,clipLimits),...            
+            getClipFrame(imp,clipLimits),...
+            getClipFrame(imfm,clipLimits),...            
             ];
     end