From 16d8ac05da3d31812f6cb2bc2a6506eb4ef477bf Mon Sep 17 00:00:00 2001
From: ac <andrew.r.cohen@drexel.edu>
Date: Sat, 31 Aug 2024 12:40:23 -0400
Subject: [PATCH] regnet corr xval qualify
---
src/MATLAB/clusterExamples/getMetaOCT.m | 31 +++++++++++
src/MATLAB/clusterExamples/goLR.m | 68 +++++++++++++++++++++++++
src/MATLAB/ncdVSvfmd/getStats.m | 19 +++++--
src/MATLAB/ncdVSvfmd/goQualify.m | 6 ++-
src/MATLAB/ncdVSvfmd/regnetCorr.m | 36 +++++++------
5 files changed, 139 insertions(+), 21 deletions(-)
create mode 100644 src/MATLAB/clusterExamples/getMetaOCT.m
create mode 100644 src/MATLAB/clusterExamples/goLR.m
diff --git a/src/MATLAB/clusterExamples/getMetaOCT.m b/src/MATLAB/clusterExamples/getMetaOCT.m
new file mode 100644
index 0000000..acf9fc5
--- /dev/null
+++ b/src/MATLAB/clusterExamples/getMetaOCT.m
@@ -0,0 +1,31 @@
+function tblMeta = getMetaOCT(ROOT,eye,scanType)
+
+tblMeta = table();
+flist = dir(fullfile(ROOT,'**/*.LEVER'));
+mx = regexp({flist.name},'P(\d+)_(.+?)_(\d+-\d+-\d+)_(\d+-\d+-\d+)_(\w\w)_(.+?)_.*.LEVER','tokens');
+tblSkip = table();
+for ff = 1:length(flist)
+ nt = table();
+
+ nt.filename = {flist(ff).name};
+ nt.folder = {flist(ff).folder};
+ if isempty(mx{ff})
+ tblSkip = [tblSkip;nt];
+ continue;
+ end
+ nt.subjectID = str2double(mx{ff}{1}{1});
+ nt.scanType = mx{ff}{1}{2};
+ nt.dx = datetime(mx{ff}{1}{3},'InputFormat','MM-dd-yyyy');
+ nt.time = {mx{ff}{1}{4}};
+ nt.eye = mx{ff}{1}{5};
+ nt.scanID = {mx{ff}{1}{6}};
+
+
+ if ~strcmp(nt.eye,eye) || ~contains(lower(scanType),lower(nt.scanType))
+ continue;
+ end
+
+ tblMeta = [tblMeta;nt];
+end
+
+tblMeta = sortrows(tblMeta,{'subjectID','filename','dx'});
\ No newline at end of file
diff --git a/src/MATLAB/clusterExamples/goLR.m b/src/MATLAB/clusterExamples/goLR.m
new file mode 100644
index 0000000..4a36f51
--- /dev/null
+++ b/src/MATLAB/clusterExamples/goLR.m
@@ -0,0 +1,68 @@
+tic
+ROOT = '/g/leverjs/Schuman_OCT/OCT/combined';
+
+PID = 'P10010';
+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'),:);
+% 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);
+
+imf = RSF.getImf(tblMeta,[[0.5,8,8] ; [5,5,5]],'');
+d = [];
+parfor i = 1:length(imf)
+ dj = [];
+ for j = 1:length(imf)
+ dj(j) = SSF.ncd_ssf_volume(imf{i},imf{j});
+ end
+ d(i,:) = dj;
+end
+
+A = Cluster.Regularize(d);
+[~,Y] = Cluster.SpectralCluster(A,3);
+
+tblMeta.Y = Y;
+
+tS = tblMeta(strcmp(tblMeta.eye,'OD'),:);
+tD = tblMeta(strcmp(tblMeta.eye,'OS'),:);
+t1 = tS;
+t2 = tD;
+
+tO = tblMeta(contains(tblMeta.scanType,'Optic'),:);
+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')
+
+plot3(t3.Y(:,1),t3.Y(:,2),t3.Y(:,3),'mx')
+plot3(t4.Y(:,1),t4.Y(:,2),t4.Y(:,3),'cs')
+legend({'ONH','macula'})
+toc
\ No newline at end of file
diff --git a/src/MATLAB/ncdVSvfmd/getStats.m b/src/MATLAB/ncdVSvfmd/getStats.m
index 3e73f24..8197425 100644
--- a/src/MATLAB/ncdVSvfmd/getStats.m
+++ b/src/MATLAB/ncdVSvfmd/getStats.m
@@ -15,12 +15,21 @@ for i = 1:length(dxx)
continue
end
+ if pData.dx(m) > pData.dx(n)
+ 4;
+ tmp = m;
+ m = n;
+ n = tmp;
+ end
if m~=n && pData.dx(m) == pData.dx(n) % same date pairs
continue
end
%
dmd = pData.md(n) - pData.md(m);
+ if dmd>0
+ continue
+ end
nt = table();
nt.subjectID = pData.subjectID(m);
idxm = find(dxx{i}(:,2)==m & dxx{i}(:,3) == m);
@@ -28,18 +37,18 @@ for i = 1:length(dxx)
nt.X = [dxx{i}(idxm,1),dxx{i}(idxn,1),dxx{i}(j,1)];
nt.velocity = (dxx{i}(j,1));
nt.md = min(pData.md(n),pData.md(m));
- nt.dmd = abs(dmd);
+ nt.dmd = dmd;
nt.mn = [m,n];
nt.ddx = pData.dx(n) - pData.dx(m);
% find sex
idx = find(tblMD.subjectID == pData.subjectID(m),1,'first');
nt.sex = tblMD.gender_full(idx);
- if nt.dmd > 4
- continue
- end
- if years(nt.ddx)>1
+ if abs(nt.dmd)>4
continue
end
+% if years(nt.ddx)>1
+% continue
+% end
% if nt.velocity > 0.9 && nt.velocity < 0.975
% continue
% end
diff --git a/src/MATLAB/ncdVSvfmd/goQualify.m b/src/MATLAB/ncdVSvfmd/goQualify.m
index f2bad0b..627ca4d 100644
--- a/src/MATLAB/ncdVSvfmd/goQualify.m
+++ b/src/MATLAB/ncdVSvfmd/goQualify.m
@@ -1,8 +1,10 @@
load('onh_kernel_corr_pt5_5_15.mat')
+% load('kcorr_VFMD_plate_blob.mat')
+% load('kernel_corr_raw.mat')
tblQualify = readtable('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.patientID = cellfun(@(x) str2double(x{1}{1}),mx,'UniformOutput',false);
for i = 1:height(pData)
idxScan = find(strcmp(tblQualify.scanID,pData.scanID{i}));
@@ -14,3 +16,5 @@ for i = 1:height(pData)
end
+tblVelo = getStats(dxx,pData,tblVFMD);
+[err,mdlX] = regnetCorr(tblVelo);
diff --git a/src/MATLAB/ncdVSvfmd/regnetCorr.m b/src/MATLAB/ncdVSvfmd/regnetCorr.m
index d7eba0c..851f350 100644
--- a/src/MATLAB/ncdVSvfmd/regnetCorr.m
+++ b/src/MATLAB/ncdVSvfmd/regnetCorr.m
@@ -1,4 +1,9 @@
-function [rmse1,rmse2] = regnetCorr(tblVelo)
+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];
t2 = tblVelo(tblVelo.dmd~=0,:);
@@ -7,22 +12,23 @@ X1 = [t2.velocity];
X2 = [t2.velocity,t2.md];
T = t2.dmd;
-% mdl1 = fitrnet(X1,T,'Activations','relu','LayerSizes',[44]);
-% cv1 = mdl1.crossval; rmse1 = cv1.kfoldLoss^0.5;
-% mdl2 = fitrnet(X2,T,'Activations','sigmoid','LayerSizes',[10]);
-% cv2 = mdl2.crossval; rmse2 = cv2.kfoldLoss^0.5;
-% 4;
-% mdl = fitrnet(X,T,"OptimizeHyperparameters","auto", "HyperparameterOptimizationOptions",struct("AcquisitionFunctionName","expected-improvement-plus"))
-
+% leave one out x-validation
err = [];
-for i = 1:length(X1)
- trainX2 = X1;
- trainX2(i) = [];
- trainY2 = T;
- trainY2(i) = [];
- mdl = fitrnet(trainX2,trainY2,'Activations','sigmoid','LayerSizes',[10]);
- pred = mdl.predict(X2(i));
+mdlX = {};
+xTarget = X1;
+for i = 1:length(xTarget)
+ trainX = xTarget;
+ trainX(i,:) = [];
+ trainY = T;
+ trainY(i) = [];
+% mdl = fitrnet(trainX2,trainY2,'Activations','sigmoid','standardize',false,'LayerSizes',[10]);
+ mdl = fitrnet(trainX,trainY,"OptimizeHyperparameters","auto", "HyperparameterOptimizationOptions",...
+ struct("AcquisitionFunctionName","expected-improvement-plus",'verbose',0,'showplots',false));
+ mdlX{i} = mdl;
+ pred = mdl.predict(xTarget(i,:));
err(i) = abs(T(i) - pred);
+ [i,mean(err)]
end
+mean(err)
4;
--
GitLab