From c1752bdda919c28b6f6e3a49a1a6bd7b37580b41 Mon Sep 17 00:00:00 2001
From: ac <andrew.r.cohen@drexel.edu>
Date: Tue, 1 Oct 2024 22:53:40 -0400
Subject: [PATCH] rsf.darktubes rnfl et al forte
---
src/MATLAB/+RSF/darkTubes.m | 25 +++++++++++++++
src/MATLAB/+RSF/getImf.m | 10 ++++--
src/MATLAB/+RSF/quantizePair.m | 6 ++--
src/MATLAB/clusterExamples/goLR.m | 33 ++++++-------------
src/MATLAB/ncdVSvfmd/getRNFL.m | 8 +++++
src/MATLAB/ncdVSvfmd/getStats.m | 29 ++++++++++++-----
src/MATLAB/ncdVSvfmd/getVFMD.m | 2 +-
src/MATLAB/ncdVSvfmd/goKernelCorr.m | 50 +++++++++--------------------
src/MATLAB/ncdVSvfmd/kernelCorr.m | 8 ++---
src/MATLAB/ncdVSvfmd/regnetCorr.m | 14 ++++----
src/MATLAB/qualify/getRNFL.m | 3 ++
src/MATLAB/qualify/goQualify.m | 33 +++++++++++++++++++
12 files changed, 138 insertions(+), 83 deletions(-)
create mode 100644 src/MATLAB/+RSF/darkTubes.m
create mode 100644 src/MATLAB/ncdVSvfmd/getRNFL.m
create mode 100644 src/MATLAB/qualify/getRNFL.m
create mode 100644 src/MATLAB/qualify/goQualify.m
diff --git a/src/MATLAB/+RSF/darkTubes.m b/src/MATLAB/+RSF/darkTubes.m
new file mode 100644
index 0000000..f65e157
--- /dev/null
+++ b/src/MATLAB/+RSF/darkTubes.m
@@ -0,0 +1,25 @@
+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;
+% vx = [1,1,1];
+sigma = 1/sqrt(3) * thickness * vx;
+
+cudaTarget = getCudaTarget;
+imf = HIP.LoG(im,sigma,cudaTarget);
+imn = imf;
+imn(imn>0) = 0;
+imn = abs(imn);
+
+bw = imbinarize(im);
+bw = bwareaopen(bw,1e3);
+bw = imdilate(bw,ones(4*thickness));
+
+imn(~bw) = 0;
+
+imfm = fibermetric(imn);
+% imfm = imfm .* imregionalmax(imfm);
+4;
+
+
diff --git a/src/MATLAB/+RSF/getImf.m b/src/MATLAB/+RSF/getImf.m
index e0cb7b7..635d6f3 100644
--- a/src/MATLAB/+RSF/getImf.m
+++ b/src/MATLAB/+RSF/getImf.m
@@ -19,13 +19,19 @@ while startTarget <= size(pData,1)
pp = spmdIndex - 1 + startTarget;
if spmdIndex <= N_CUDA_PROCS && pp <= size(pData,1)
strDB = fullfile(pData.folder{pp},pData.filename{pp});
- [im] = SSF.loadImage(strDB,1);
+ % RSF images and metadata via lever file
+ [im,CONSTANTS] = leversc.loadImage(strDB,1,1);
% 3-channel e.g. plate,blob,0
for c = 1:length(radii)
if all(0==radii{c})
im1(:,:,:,c) = zeros(size(im));
else
- im1(:,:,:,c) = RSF.imPreProcess(im,radii{c});
+ if length(radii{c}) == 1
+ % hardcode for darktubes...careful...ACK
+ im1(:,:,:,c) = RSF.darkTubes(im,CONSTANTS,radii{c});
+ else
+ im1(:,:,:,c) = RSF.imPreProcess(im,radii{c});
+ end
end
end
end
diff --git a/src/MATLAB/+RSF/quantizePair.m b/src/MATLAB/+RSF/quantizePair.m
index 7bda641..03dc8e5 100644
--- a/src/MATLAB/+RSF/quantizePair.m
+++ b/src/MATLAB/+RSF/quantizePair.m
@@ -1,9 +1,9 @@
% clipLimits = [minVal,maxVal];
function [i1q,i2q] = quantizePair(i1,i2)
-if isa(i1,'uint8') && isa(i2,'uint8')
- return
-end
+% if isa(i1,'uint8') && isa(i2,'uint8')
+% return
+% end
pix = [i1(find(i1));i2(find(i2))];
diff --git a/src/MATLAB/clusterExamples/goLR.m b/src/MATLAB/clusterExamples/goLR.m
index 4c47c8b..8fa05c2 100644
--- a/src/MATLAB/clusterExamples/goLR.m
+++ b/src/MATLAB/clusterExamples/goLR.m
@@ -1,5 +1,5 @@
tic
-ROOT = '/g/leverjs/Schuman_OCT/OCT/combined';
+ROOT = '/g/leverjs/Schuman_OCT/OCT/qualified';
PID = 'P10010';
flist = dir(fullfile(ROOT,['/' PID '*.LEVER']));
@@ -13,25 +13,12 @@ tblMeta = tblMeta(contains(tblMeta.scanType,'Optic'),:);
% tblMeta = tblMeta(contains(tblMeta.scanType,'Macula'),:);
% tblMeta = tblMeta(contains(tblMeta.eye,'OD'),:);
-tblQualify = readtable('../ncdVSvfmd/nyu onh matched.xlsx');
-mx = regexp(tblQualify.file_name,'P(\d+)_(.+?)_(\d+-\d+-\d+)_(\d+-\d+-\d+)_(\w\w)_(.+?)_.*.img','tokens');
-tblQualify.scanID = cellfun(@(x) x{1}{6},mx,'UniformOutput',false);
-tblQualify.PID = cellfun(@(x) str2double(x{1}{1}),mx,'UniformOutput',false);
-PID = str2double(PID(2:end));
-idx = find(cell2mat(tblQualify.PID) == PID);
-
-for i = 1:height(tblMeta)
- idx = find(strcmp(tblQualify.scanID,tblMeta.scanID{i}));
- if any(strcmp(tblQualify.qualification_status(idx),'No'))
- tblMeta.bDisqualified(i) = true;
- end
-end
-tblMeta(tblMeta.bDisqualified,:) = [];
p = ljsStartParallel(96);
-radii = {[0.5,8,8],[5,5,5],[0,0,0]}; % plate,blob
+% radii = {[0.5,8,8],[5,5,5],[0,0,0]}; % plate,blob
% radii = {[0.5,8,8]}; % plate
+radii = {[3]} % dark tubes
imf = RSF.getImf(tblMeta,radii);
d = [];
parfor i = 1:length(imf)
@@ -57,13 +44,13 @@ tM = tblMeta(contains(tblMeta.scanType,'Macular'),:);
t3 = tO;
t4 = tM;
-clf;hold on
-plot3(t1.Y(:,1),t1.Y(:,2),t1.Y(:,3),'r*')
-plot3(t2.Y(:,1),t2.Y(:,2),t2.Y(:,3),'og')
-legend({'OS','OD'})
-xlabel('NCD1')
-ylabel('NCD2')
-zlabel('NCD3')
+ clf;hold on
+ plot3(t1.Y(:,1),t1.Y(:,2),t1.Y(:,3),'r*')
+ plot3(t2.Y(:,1),t2.Y(:,2),t2.Y(:,3),'og')
+ legend({'OS','OD'})
+ xlabel('NCD1')
+ ylabel('NCD2')
+ zlabel('NCD3')
plot3(t3.Y(:,1),t3.Y(:,2),t3.Y(:,3),'mx')
plot3(t4.Y(:,1),t4.Y(:,2),t4.Y(:,3),'cs')
diff --git a/src/MATLAB/ncdVSvfmd/getRNFL.m b/src/MATLAB/ncdVSvfmd/getRNFL.m
new file mode 100644
index 0000000..fb7bb82
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/getRNFL.m
@@ -0,0 +1,8 @@
+function iRNFL = getRNFL(p1,tblRNFL)
+
+% find p1 in tblRNFL and return avg thickness
+[~,f1,~] = fileparts(p1.filename);
+[~,tblFname,~] = fileparts(tblRNFL.file_name);
+idx1 = find(strcmp(f1,tblFname),1,'first');
+iRNFL = tblRNFL{idx1,"AVERAGETHICKNESS"};
+4;
diff --git a/src/MATLAB/ncdVSvfmd/getStats.m b/src/MATLAB/ncdVSvfmd/getStats.m
index 8197425..9a15882 100644
--- a/src/MATLAB/ncdVSvfmd/getStats.m
+++ b/src/MATLAB/ncdVSvfmd/getStats.m
@@ -1,19 +1,21 @@
-function tblVelo = getStats(dxx, pData, tblMD)
-if ~exist('idxDel','var')
- idxDel = [];
+function tblVelo = getStats(dxx, pData, tblMD, tblRNFL)
+
+if ~exist('tblRNFL','var')
+ tblRNFL = [];
end
+
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
% check first if either is disqualified
- if pData.disqualified(m) || pData.disqualified(n)
- continue
- end
+
if pData.dx(m) > pData.dx(n)
4;
@@ -30,26 +32,37 @@ for i = 1:length(dxx)
if dmd>0
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);
idxn = find(dxx{i}(:,2)==n & dxx{i}(:,3) == n);
nt.X = [dxx{i}(idxm,1),dxx{i}(idxn,1),dxx{i}(j,1)];
- nt.velocity = (dxx{i}(j,1));
+ nt.ncd = (dxx{i}(j,1));
nt.md = min(pData.md(n),pData.md(m));
nt.dmd = dmd;
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);
if abs(nt.dmd)>4
continue
end
+
+% on the unqualified data vs. qualified data
% if years(nt.ddx)>1
% continue
% end
-% if nt.velocity > 0.9 && nt.velocity < 0.975
+% if nt.ncd > 0.9 && nt.ncd< 0.975
% continue
% end
if m ~= n
diff --git a/src/MATLAB/ncdVSvfmd/getVFMD.m b/src/MATLAB/ncdVSvfmd/getVFMD.m
index 92e08e5..4626b63 100644
--- a/src/MATLAB/ncdVSvfmd/getVFMD.m
+++ b/src/MATLAB/ncdVSvfmd/getVFMD.m
@@ -10,7 +10,7 @@ fx = {'nyu VF qualified data great than -6 data export.xlsx',...
tblVFMD = table();
for ff = 1 : length(fx)
- tx = readtable(fullfile("importMetadata/dataExcel/",fx{ff}));
+ tx = readtable(fullfile("../importMetadata/dataExcel/",fx{ff}));
% filter 1
idxVFMD = contains(lower(tx.qualification_status),'yes');
tx = tx(idxVFMD,:);
diff --git a/src/MATLAB/ncdVSvfmd/goKernelCorr.m b/src/MATLAB/ncdVSvfmd/goKernelCorr.m
index 992b050..9451fd9 100644
--- a/src/MATLAB/ncdVSvfmd/goKernelCorr.m
+++ b/src/MATLAB/ncdVSvfmd/goKernelCorr.m
@@ -1,42 +1,22 @@
datetime
kernelCorrStart = tic();
eye = 'OD';
-ROOT = '/g/leverjs/Schuman_OCT/OCT/combined';
+ROOT = '/g/leverjs/Schuman_OCT/OCT/qualified';
% scanType = 'Macular Cube 200x200';
scanType = 'Optic Disc Cube 200x200';
load('vfmd_meta.mat');
-[rho2,tblVelo, pData, dxx] = kernelCorr([0.5,8,8],tblVFMD,tblMeta,scanType,eye);
-[rmse1,rmse2] = regnetCorr(tblVelo)
-%
-% tblResults = table();
-% for rx = 1:length(radii)
-% rr = radii{rx};
-% [rho2,tblVelo, pData, dxx] = kernelCorr(rr,tblVFMD,tblMeta,scanType,eye);
-% nt = table();
-% nt.rr = {rr};
-% nt.tblVelo = {tblVelo};
-% nt.pData = {pData};
-% nt.dxx = {dxx};
-% nt.r2 = rho2;
-% if iscell(rr)
-% rstr = cellfun(@mat2str,rr,'UniformOutput',false);
-% rstr = horzcat(rstr{:});
-% else
-% rstr = mat2str(rr);
-% end
-% nt.rstr = {rstr};
-% tblResults = [tblResults;nt]
-% end
-%
-% for i = 1:height(tblResults)
-% [rmse1,rmse2] = regnetCorr(tblResults.tblVelo{i});
-% tblResults.rmse(i,:) = [rmse1,rmse2];
-% tblVelo = tblResults.tblVelo{i};
-% x = tblVelo.velocity;
-% y = tblVelo.dmd;
-% yt = log10(1+y);
-% [f,gofX]=fit(x,yt,'poly1');
-% [f,gof]=fit(x,y,'poly1');
-% tblResults.r2a(i,1:2) = [gof.adjrsquare,gofX.adjrsquare]
-% end
+% tblVFMD = getVFMD(eye);
+tblMeta = getMetaOCT(ROOT,eye,scanType);
+
+
+load('../qualify/onh_rnfl_thickness.mat');
+
+radii = {[0.5,8,8],[5,5,5],[4]}; % plate,blob
+% radii = {[0.5,8,8]}; % plate
+
+
+[rho2,tblVelo, pData, dxx] = kernelCorr(radii,tblVFMD,tblMeta,tblRNFL);
+
+
+[rmse1,~] = regnetCorr(tblVelo);median(rmse1)
diff --git a/src/MATLAB/ncdVSvfmd/kernelCorr.m b/src/MATLAB/ncdVSvfmd/kernelCorr.m
index 82dbcc2..1d313ce 100644
--- a/src/MATLAB/ncdVSvfmd/kernelCorr.m
+++ b/src/MATLAB/ncdVSvfmd/kernelCorr.m
@@ -1,4 +1,4 @@
-function [rho2,tblVelo, pData, dxx] = kernelCorr(radii,tblVFMD,tblMeta,scanType,eye)
+function [rho2,tblVelo, pData, dxx] = kernelCorr(radii,tblVFMD,tblMeta,tblRNFL)
kernelCorrStart = tic();
@@ -40,7 +40,7 @@ end
p=ljsStartParallel(96);
-imf = RSF.getImf(pData,radii,scanType);
+imf = RSF.getImf(pData,radii);
[C,ia,ic] = unique(pData.subjectID);
@@ -84,8 +84,8 @@ parfor wx = 1:length(workList)
end
end
-tblVelo = getStats(dxx,pData,tblVFMD);
-[rho,pCorr] = corr((tblVelo.velocity),(tblVelo.dmd));
+tblVelo = getStats(dxx,pData,tblVFMD,tblRNFL);
+[rho,pCorr] = corr((tblVelo.ncd),(tblVelo.dmd));
rho2 = rho^2
toc(kernelCorrStart)
4;
diff --git a/src/MATLAB/ncdVSvfmd/regnetCorr.m b/src/MATLAB/ncdVSvfmd/regnetCorr.m
index 851f350..6a01931 100644
--- a/src/MATLAB/ncdVSvfmd/regnetCorr.m
+++ b/src/MATLAB/ncdVSvfmd/regnetCorr.m
@@ -1,21 +1,21 @@
function [err,mdlX] = regnetCorr(tblVelo)
-% tblVelo.velocity is NCD
% tblVelo.md is vfmd for first of two image pairs
% tblVelo.dmd is delta VFMD between two image pairs
% tblVelo = getStats(dxx,pData,[],tblVFMD);
-% X = [tblVelo.velocity,tblVelo.md];
+% X = [tblVelo.ncd,tblVelo.md];
t2 = tblVelo(tblVelo.dmd~=0,:);
-X1 = [t2.velocity];
-X2 = [t2.velocity,t2.md];
+X1 = [t2.ncd];
+X2 = [t2.ncd,t2.md];
T = t2.dmd;
+% T = t2.deltaRNFL;
% leave one out x-validation
err = [];
mdlX = {};
-xTarget = X1;
+xTarget = X2;
for i = 1:length(xTarget)
trainX = xTarget;
trainX(i,:) = [];
@@ -27,8 +27,8 @@ for i = 1:length(xTarget)
mdlX{i} = mdl;
pred = mdl.predict(xTarget(i,:));
err(i) = abs(T(i) - pred);
- [i,mean(err)]
+ [i,median(err)]
end
-mean(err)
+% mean(err)
4;
diff --git a/src/MATLAB/qualify/getRNFL.m b/src/MATLAB/qualify/getRNFL.m
new file mode 100644
index 0000000..50f0dc4
--- /dev/null
+++ b/src/MATLAB/qualify/getRNFL.m
@@ -0,0 +1,3 @@
+% getRNFL
+tx_onh = readtable('nyu onh matched.xlsx');
+save('onh_rnfl_thickness.mat','tx_onh');
diff --git a/src/MATLAB/qualify/goQualify.m b/src/MATLAB/qualify/goQualify.m
new file mode 100644
index 0000000..1c94df3
--- /dev/null
+++ b/src/MATLAB/qualify/goQualify.m
@@ -0,0 +1,33 @@
+ROOT = '/g/leverjs/Schuman_OCT/OCT/combined';
+target = '/g/leverjs/Schuman_OCT/OCT/qualified';
+if ~exist(target,'dir')
+ mkdir(target);
+end
+% flist = dir(fullfile(ROOT,'**/*.LEVER'));
+% qONH = readtable('nyu onh matched.xlsx');
+% qMacula = readtable('nyu mac matched.xlsx');
+% [~,fONH,~] = fileparts(qONH.file_name);
+% [~,fMacula,~] = fileparts(qMacula.file_name);
+qlist = [];
+for i = 1 : length(flist)
+
+ [~,fname,~] = fileparts(flist(i).name);
+ idx = find(strcmp(fname,fMacula));
+ if ~isempty(idx)
+ tidx = qMacula(idx,:);
+ else
+ idx = find(strcmp(fname,fONH));
+ tidx = qONH(idx,:);
+ end
+ if isempty(idx)
+ continue
+ end
+ qs = lower(tidx.qualification_status);
+ if any(strcmp(qs,'no'))
+ continue
+ end
+ qlist = [qlist,flist(i)];
+ src = fullfile(flist(i).folder,fname);
+ cmd = ['cp "' src '".* ' target];
+ system(cmd);
+end
--
GitLab