diff --git a/src/MATLAB/+RSF/getImf.asv b/src/MATLAB/+RSF/getImf.asv
new file mode 100644
index 0000000000000000000000000000000000000000..e7d53f73d153bacc5b528492b5f0f2e156799eed
--- /dev/null
+++ b/src/MATLAB/+RSF/getImf.asv
@@ -0,0 +1,86 @@
+function imf = getImf(pData,radii,scanType)
+
+% radii is 1-element or 3-element cell array. [0,0,0] for unused channel
+% FLIF takes 1 or 3 UBYTE
+radii = {[0.5,8,8],[5,5,5],[10,10,10]}; % plate,blob
+
+% NGPU = HIP.Cuda.DeviceCount;
+% p = ljsStartParallel(8*NGPU);
+
+N_CUDA_PROCS = 48;
+startTarget = 1;
+imf = {};
+im1 = Composite();
+pp = Composite();
+while startTarget <= size(pData,1)
+ spmd
+ 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);
+ % 3-channel 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});
+ end
+ end
+ end
+ startTarget = startTarget + N_CUDA_PROCS;
+ for i = 1:N_CUDA_PROCS
+ if pp{i} <= size(pData,1)
+ imf{pp{i}} = im1{i};
+ end
+ end
+end
+
+px = {}; % critical voxels
+mx = []; % mean by stack
+sx = []; % sigma by stack
+for i = 1 : length(imf)
+ for c = 1:size(imf{i},4)
+ i1 = imf{i}(:,:,:,c);
+ px{i,c} = i1(find(i1));
+ mx(i,c) = mean(px{i,c});
+ sx(i,c) = std(px{i,c});
+ end
+end
+
+% clipRange = [2.5,97.5];
+clipRange = [0.5,99.5];
+for c = 1:size(imf{1},4)
+ p1 = vertcat(px{:,c});
+ px_pos = p1(p1>0);
+ px_neg = abs(p1(p1<0));
+
+ clip_pos(c,:) = prctile(px_pos,clipRange);
+ clip_neg(c,:) = prctile(px_neg,clipRange);
+end
+4;
+for i = 1:length(imf)
+ for c = 1 : size(imf{i},4)
+ if all(0==imf{i}(:,:,c))
+ continue
+ end
+ imp = imf{i}(:,:,:,c);
+ imp(imp<0)=0;
+ imn = imf{i}(:,:,:,c);
+ imn(imn>0)=0;
+ imn=abs(imn);
+ imq_pos = SSF.quantize8(imp,clip_pos(c,:));
+ imq_neg = SSF.quantize8(imn,clip_neg(c,:));
+ if all(0==imq_pos)
+ imf{i}(:,:,:,c) = imq_neg;
+ else
+ % [0,127] for imq_neg, [128,255] imq_pos
+ imq_pos = (0.5 .* imq_pos) + uint8(128 .* (imq_pos>0));
+ imq_neg = 0.5 .* imq_neg ;
+ imf{i}(:,:,:,c) = (imq_neg + imq_pos);
+ end
+
+ end
+
+ imf{i} = uint8(imf{i});
+end
+4;
diff --git a/src/MATLAB/+RSF/getImf.m b/src/MATLAB/+RSF/getImf.m
index b77b44e2f5fe0f17aa0161b9c470554893d6df31..67e2384714a6eeb4916ea20850a4b7e66b8d4aa0 100644
--- a/src/MATLAB/+RSF/getImf.m
+++ b/src/MATLAB/+RSF/getImf.m
@@ -1,45 +1,88 @@
function imf = getImf(pData,radii,scanType)
-NGPU = HIP.Cuda.DeviceCount;
-p = ljsStartParallel(8*NGPU);
+% radii is 1-element or 3-element cell array. [0,0,0] for unused channel
+% FLIF takes 1 or 3 UBYTE
+radii = {[0.5,8,8],[5,5,5],[10,10,10]}; % plate,blob
+% radii = {[0.5,8,8]}; % plate
+% NGPU = HIP.Cuda.DeviceCount;
+% p = ljsStartParallel(8*NGPU);
+
+N_CUDA_PROCS = 48;
+startTarget = 1;
imf = {};
-px = {};
-parfor pp = 1:height(pData)
- strDB = fullfile(pData.folder{pp},pData.filename{pp});
- [im] = SSF.loadImage(strDB,1);
- imf{pp} = RSF.imPreProcess(im,radii);
- p1 = imf{pp}(:);
- p1 = p1(find(p1));
- px{pp} = p1;
+im1 = Composite();
+pp = Composite();
+while startTarget <= size(pData,1)
+ spmd
+ 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);
+ % 3-channel 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});
+ end
+ end
+ end
+ end
+ startTarget = startTarget + N_CUDA_PROCS;
+ for i = 1:N_CUDA_PROCS
+ if pp{i} <= size(pData,1)
+ imf{pp{i}} = im1{i};
+ end
+ end
end
-px = px{:};
-px_pos = px(px>0);
-px_neg = abs(px(px<0));
+% px = {}; % critical voxels
+% mx = []; % mean by stack
+% sx = []; % sigma by stack
+% for i = 1 : length(imf)
+% for c = 1:size(imf{i},4)
+% i1 = imf{i}(:,:,:,c);
+% px{i,c} = i1(find(i1));
+% mx(i,c) = mean(px{i,c});
+% sx(i,c) = std(px{i,c});
+% end
+% end
-% clip_pos = [mean(px_pos)-std(px_pos), mean(px_pos)+std(px_pos)];
-% clip_neg = [mean(px_neg)-std(px_neg), mean(px_neg)+std(px_neg)];
-% clip_pos = [min(px_pos),max(px_pos)];
-% clip_neg = [min(px_neg),max(px_neg)];
-clip_pos = prctile(px_pos,[0.5,99.5]);
-clip_neg = prctile(px_neg,[0.5,99.5]);
-
-for i = 1:length(imf)
- imp = imf{i};
- imp(imp<0)=0;
- imn = imf{i};
- imn(imn>0)=0;
- imn=abs(imn);
- imq_pos = SSF.quantize8(imp,clip_pos);
- imq_neg = SSF.quantize8(imn,clip_neg);
- if all(0==imq_pos)
- imf{i} = imq_neg;
- else
- % [0,127] for imq_neg, [128,255] imq_pos
- imq_pos = (0.5 .* imq_pos) + uint8(128 .* (imq_pos>0));
- imq_neg = 0.5 .* imq_neg ;
- imf{i} = imq_neg + imq_pos;
- end
-end
-4;
+% clipRange = [2.5,97.5];
+% % clipRange = [0.5,99.5];
+% for c = 1:size(imf{1},4)
+% p1 = vertcat(px{:,c});
+% px_pos = p1(p1>0);
+% px_neg = abs(p1(p1<0));
+%
+% clip_pos(c,:) = prctile(px_pos,clipRange);
+% clip_neg(c,:) = prctile(px_neg,clipRange);
+% end
+% 4;
+% for i = 1:length(imf)
+% for c = 1 : size(imf{i},4)
+% if all(0==imf{i}(:,:,c))
+% continue
+% end
+% imp = imf{i}(:,:,:,c);
+% imp(imp<0)=0;
+% imn = imf{i}(:,:,:,c);
+% imn(imn>0)=0;
+% imn=abs(imn);
+% imq_pos = SSF.quantize8(imp,clip_pos(c,:));
+% imq_neg = SSF.quantize8(imn,clip_neg(c,:));
+% if all(0==imq_pos)
+% imf{i}(:,:,:,c) = imq_neg;
+% else
+% % [0,127] for imq_neg, [128,255] imq_pos
+% imq_pos = (0.5 .* imq_pos) + uint8(128 .* (imq_pos>0));
+% imq_neg = 0.5 .* imq_neg ;
+% imf{i}(:,:,:,c) = (imq_neg + imq_pos);
+% end
+%
+% end
+%
+% imf{i} = uint8(imf{i});
+% end
+% 4;
diff --git a/src/MATLAB/+RSF/quantize8.asv b/src/MATLAB/+RSF/quantize8.asv
new file mode 100644
index 0000000000000000000000000000000000000000..f7c98e1698f78f837f21417ab36d473c59d59331
--- /dev/null
+++ b/src/MATLAB/+RSF/quantize8.asv
@@ -0,0 +1,18 @@
+% clipLimits = [minVal,maxVal];
+function [im, clipLimits] = quantize8(i1,i2)
+
+if isa(i1,'uint8') && isa(i2,'uint8')
+ return
+end
+
+pix = [i1(find(i1)),i2(find(i2))];
+% clip at 99.9%
+clipLimits = prctile(pix,[0.05,99.95]);
+
+im = max(im,clipLimits(1));
+im = min(im,clipLimits(end));
+im = (im - clipLimits(1)) ./ (clipLimits(end) - clipLimits(1));
+im = im2uint8(im);
+im(0 == im_in) = 0;
+4;
+
diff --git a/src/MATLAB/+RSF/quantizePair.m b/src/MATLAB/+RSF/quantizePair.m
new file mode 100644
index 0000000000000000000000000000000000000000..7bda6410bb87bbc71dd82d2e3a096cb961337db6
--- /dev/null
+++ b/src/MATLAB/+RSF/quantizePair.m
@@ -0,0 +1,26 @@
+% clipLimits = [minVal,maxVal];
+function [i1q,i2q] = quantizePair(i1,i2)
+
+if isa(i1,'uint8') && isa(i2,'uint8')
+ return
+end
+
+pix = [i1(find(i1));i2(find(i2))];
+
+% clip at 99.9%
+clipLimits = prctile(pix,[0.05,99.95]);
+
+i1q = applyClip(i1,clipLimits);
+i2q = applyClip(i2,clipLimits);
+
+4;
+
+function im = applyClip(im,clipLimits)
+im_in = im; % remember the zeros!
+im = max(im,clipLimits(1));
+im = min(im,clipLimits(end));
+im = (im - clipLimits(1)) ./ (clipLimits(end) - clipLimits(1));
+im = uint8(254 * im + 1); % 8 bit quantization [0,1]->[1,255]
+im(0 == im_in) = 0;
+4;
+
diff --git a/src/MATLAB/metadata/goLR.m b/src/MATLAB/metadata/goLR.m
index 3879125c9b184def8b8eb62581f48483f6a806ff..5672154465203e5a4f4caef7ba96c838ffd03840 100644
--- a/src/MATLAB/metadata/goLR.m
+++ b/src/MATLAB/metadata/goLR.m
@@ -1,16 +1,6 @@
tic
ROOT = '/g/leverjs/Schuman_OCT/OCT/03-22-2024';
-% flist = dir(fullfile(ROOT,'**/*.LEVER'));
-% id = regexp({flist.name},'P(\d)+','tokens');
-% id = cellfun(@(x) str2double(x{1}),id);
-% [uid,ia,ic] = unique(id);
-%
-% hx = [];
-% for i = 1:length(uid)
-% hx(i) = length(find(id==uid(i)));
-% end
-
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');
@@ -19,15 +9,12 @@ tblMeta.PID = str2double(tblMeta.PID);
tblMeta.folder = {flist.folder}';
tblMeta.filename = {flist.name}';
-tblMeta = tblMeta(contains(tblMeta.scanType,'Optic'),:);
+% tblMeta = tblMeta(contains(tblMeta.scanType,'Optic'),:);
+tblMeta = tblMeta(contains(tblMeta.scanType,'Macula'),:);
% tblMeta = tblMeta(contains(tblMeta.eye,'OD'),:);
-imf = RSF.getImf(tblMeta,[0.5,8,8],'');
-imf2 = RSF.getImf(tblMeta,[5,5,5],'');
-for i = 1:length(imf)
- imf{i}(:,:,:,2) = imf2{i};
- imf{i}(:,:,:,3) = 0 .* imf2{i};
-end
-p = ljsStartParallel(64);
+p = ljsStartParallel(96);
+
+imf = RSF.getImf(tblMeta,[[0.5,8,8] ; [5,5,5]],'');
d = [];
parfor i = 1:length(imf)
dj = [];
@@ -47,15 +34,16 @@ 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;
+% 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')
-plot3(t3.Y(:,1),t3.Y(:,2),t3.Y(:,3),'mx')
-plot3(t4.Y(:,1),t4.Y(:,2),t4.Y(:,3),'cs')
+
+% plot3(t3.Y(:,1),t3.Y(:,2),t3.Y(:,3),'mx')
+% plot3(t4.Y(:,1),t4.Y(:,2),t4.Y(:,3),'cs')
toc
\ No newline at end of file
diff --git a/src/MATLAB/metadata/leftVright_3_3_3_d.mat b/src/MATLAB/metadata/leftVright_3_3_3_d.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e8c2b5dcc20f3bd5c2bbad88466a6fff5b2f11f0
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_3_3_3_d.mat differ
diff --git a/src/MATLAB/metadata/leftVright_5_5_5_d.mat b/src/MATLAB/metadata/leftVright_5_5_5_d.mat
new file mode 100644
index 0000000000000000000000000000000000000000..015fa8b8a1368c302731d40fb34fd36be5149c36
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_5_5_5_d.mat differ
diff --git a/src/MATLAB/metadata/leftVright_8_8_8_d.mat b/src/MATLAB/metadata/leftVright_8_8_8_d.mat
new file mode 100644
index 0000000000000000000000000000000000000000..57742ab1a6ee57568aa8f92a4a2c6b56c97dfde2
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_8_8_8_d.mat differ
diff --git a/src/MATLAB/metadata/leftVright_MAC_ONH_plate_blob.mat b/src/MATLAB/metadata/leftVright_MAC_ONH_plate_blob.mat
new file mode 100644
index 0000000000000000000000000000000000000000..244895e9cd0412e37b85d92b4cd48b4334341e6f
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_MAC_ONH_plate_blob.mat differ
diff --git a/src/MATLAB/metadata/leftVright_OD_OS_plate_blob.mat b/src/MATLAB/metadata/leftVright_OD_OS_plate_blob.mat
new file mode 100644
index 0000000000000000000000000000000000000000..945f0f87208afcfcf7c35011e92f8721df52035e
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_OD_OS_plate_blob.mat differ
diff --git a/src/MATLAB/metadata/leftVright_OD_OS_plate_blob_macula.mat b/src/MATLAB/metadata/leftVright_OD_OS_plate_blob_macula.mat
new file mode 100644
index 0000000000000000000000000000000000000000..11aef006c4b53ee4302b9c5506b8795f61c7dec7
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_OD_OS_plate_blob_macula.mat differ
diff --git a/src/MATLAB/metadata/leftVright_pt5_10_10_d.mat b/src/MATLAB/metadata/leftVright_pt5_10_10_d.mat
new file mode 100644
index 0000000000000000000000000000000000000000..4b2f570003e28261a4524e026bcbf788a3b84e33
Binary files /dev/null and b/src/MATLAB/metadata/leftVright_pt5_10_10_d.mat differ
diff --git a/src/MATLAB/ncdVSvfmd/getMetaOCT.m b/src/MATLAB/ncdVSvfmd/getMetaOCT.m
new file mode 100644
index 0000000000000000000000000000000000000000..acf9fc5352cd522b236a7feea224c9c1a2d8d672
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/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/ncdVSvfmd/getStats.m b/src/MATLAB/ncdVSvfmd/getStats.m
new file mode 100644
index 0000000000000000000000000000000000000000..56b12eb542b7ecab65d31adc089757bd54c9397e
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/getStats.m
@@ -0,0 +1,54 @@
+
+function tblVelo = getStats(dxx, pData, idxDel, tblMD)
+if ~exist('idxDel','var')
+ idxDel = [];
+end
+tblVelo = table();
+idxUnused = 1:height(pData);
+for i = 1:length(dxx)
+ for j = 1:size(dxx{i},1)
+ m = dxx{i}(j,2);
+ n = dxx{i}(j,3);
+ if ~isempty(intersect(m,idxDel)) || ~isempty(intersect(n,idxDel))
+ continue
+ end
+
+
+ if m~=n && pData.dx(m) == pData.dx(n) % same date pairs
+ continue
+ end
+%
+
+ dmd = pData.md(n) - pData.md(m);
+ 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.md = min(pData.md(n),pData.md(m));
+ nt.dmd = abs(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
+ continue
+ end
+% if nt.velocity > 0.9 && nt.velocity < 0.975
+% continue
+% end
+ if m ~= n
+ idxUnused = setdiff(idxUnused,m);
+ idxUnused = setdiff(idxUnused,n);
+ end
+ tblVelo = [tblVelo;nt];
+ end
+end
+
+% tblVelo = trimSelfSimilar(tblVelo);
+4;
diff --git a/src/MATLAB/ncdVSvfmd/getVFMD.m b/src/MATLAB/ncdVSvfmd/getVFMD.m
new file mode 100644
index 0000000000000000000000000000000000000000..92e08e5386fdad126d4713b353c322c7cf318dec
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/getVFMD.m
@@ -0,0 +1,67 @@
+function tblVFMD = getVFMD(eye)
+
+if ~exist('eye','var')
+ eye = 'OD';
+end
+fx = {'nyu VF qualified data great than -6 data export.xlsx',...
+ 'nyu VF qualified data from -12 to -6 data export.xlsx',...
+ 'nyu VF qualified data less than -12 data export.xlsx'};
+% read from 3 files
+tblVFMD = table();
+for ff = 1 : length(fx)
+
+ tx = readtable(fullfile("importMetadata/dataExcel/",fx{ff}));
+ % filter 1
+ idxVFMD = contains(lower(tx.qualification_status),'yes');
+ tx = tx(idxVFMD,:);
+ % filter 2
+ idxVFMD = contains(lower(tx.threshold),lower('Visual Field 24-2 Test Pattern'));
+ tx = tx(idxVFMD,:);
+ % filter 3
+ for i = 1:height(tx)
+ tx.operatorError(i) = sum([tx.false_positive_p(i), tx.false_negative_p(i), tx.fixation_loss_c_p(i)]);
+ end
+ idxVFMD = tx.operatorError < 33;
+ tx = tx(idxVFMD,:);
+
+ % filter by eye
+ idxVFMD = contains(lower(tx.Eye),lower(eye));
+ tx = tx(idxVFMD,:);
+
+ tblVFMD = [tblVFMD;tx];
+end
+
+tblVFMD.subjectID = cellfun(@(x) str2double(x),tblVFMD.project_id);
+tblVFMD.dx = tblVFMD.visit_date_date;
+tblVFMD = sortrows(tblVFMD,{'subjectID','dx'});
+tblVFMD = unique(tblVFMD,'rows');
+
+% filter by date ,patient
+[pC, p_ia, p_ic] = unique(tblVFMD.subjectID);
+for pp = 1:length(pC)
+ idxVFMD = tblVFMD.subjectID == pC(pp);
+ tbl1 = tblVFMD(idxVFMD,:);
+ dx = unique(tbl1.visit_date_date);
+ idxDelAll = [];
+ for dd = 1:length(dx)
+ d1 = tbl1(tbl1.visit_date_date == dx(dd),:);
+ if 1 == height(d1)
+ continue
+ end
+ % need to choose best visit
+ px = [];
+ for vv = 1:height(d1)
+ px(vv,:) = [d1.false_positive_p(vv),d1.false_negative_p(vv),d1.fixation_loss_c_p(vv)];
+ end
+ px = sum(abs(px),2);
+ [~,idxMin] = min(px);
+ idxDel = find(idxVFMD & tblVFMD.dx == dx(dd));
+ idxDel(idxMin) = []; % preserve the most likely
+ idxDelAll = [idxDelAll,idxDel];
+
+ 4;
+ end
+ tblVFMD(idxDelAll,:) = []; % remove the least likely
+end
+tblVFMD.mdString = tblVFMD.md;
+tblVFMD.md = cellfun(@str2double,tblVFMD.md);
diff --git a/src/MATLAB/ncdVSvfmd/goKernelCorr.m b/src/MATLAB/ncdVSvfmd/goKernelCorr.m
new file mode 100644
index 0000000000000000000000000000000000000000..992b050f2472db01dd10eaef12a4cf0dc0e3c587
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/goKernelCorr.m
@@ -0,0 +1,42 @@
+datetime
+kernelCorrStart = tic();
+eye = 'OD';
+ROOT = '/g/leverjs/Schuman_OCT/OCT/combined';
+% 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
diff --git a/src/MATLAB/ncdVSvfmd/kernelCorr.m b/src/MATLAB/ncdVSvfmd/kernelCorr.m
new file mode 100644
index 0000000000000000000000000000000000000000..8320c86c516c154329814d3781f560a91e0f8978
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/kernelCorr.m
@@ -0,0 +1,93 @@
+function [rho2,tblVelo, pData, dxx] = kernelCorr(radii,tblVFMD,tblMeta,scanType,eye)
+
+kernelCorrStart = tic();
+
+mdInfo = [];
+pid = unique(tblVFMD.subjectID);
+for i = 1:length(pid)
+ idx = find(tblVFMD.subjectID == pid(i));
+ visits = unique(tblVFMD.dx(idx));
+ if length(visits)<2,continue,end
+ md1 = [];
+ md1.id = pid(i);
+ md1.dx = tblVFMD.dx(idx);
+ md1.md = tblVFMD.md(idx);
+ mdInfo = [mdInfo,md1];
+end
+
+for i = 1:length(mdInfo)
+ idx = find(tblMeta.subjectID == mdInfo(i).id);
+ [C,ia,ib] = intersect(mdInfo(i).dx,tblMeta.dx(idx));
+ mdInfo(i).dxOCT = C;
+ mdInfo(i).dx = mdInfo(i).dx(ia);
+ mdInfo(i).md = mdInfo(i).md(ia);
+end
+
+octLen = arrayfun(@(x) length(x.dxOCT),mdInfo);
+mdInfo(octLen<2) = [];
+
+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.md = repmat(mdInfo(i).md(pj),height(p1),1);
+ pData = [pData;p1];
+ end
+end
+
+% pData = pData(1:20,:); % shorten for test
+
+p=ljsStartParallel(96);
+
+imf = RSF.getImf(pData,radii,scanType);
+
+[C,ia,ic] = unique(pData.subjectID);
+
+
+workList = [];
+for px = 1:length(C)
+ work1 = [];
+ idx = find(ic==px);
+
+ work1.imf = imf(idx);
+
+ wl = [idx,idx];
+ wl = [wl;nchoosek(idx,2)];
+ work1.wl = wl;
+
+ idx0 = [1:length(idx)]';
+ wl0 = [idx0,idx0];
+ wl0 = [wl0;nchoosek(idx0,2)];
+ work1.wl0 = wl0;
+
+ workList = [workList,work1];
+end
+
+
+dxx = {};
+parfor wx = 1:length(workList)
+ for ii = 1:size(workList(wx).wl0,1)
+ m = workList(wx).wl0(ii,1) ;
+ n = workList(wx).wl0(ii,2);
+ im_n = workList(wx).imf{n};
+ im_m = workList(wx).imf{m};
+ d1 = NCD.ncd_ssf_volume(im_m,im_n);
+ if m==n
+ d2 = d1;
+ else
+ d2 = NCD.ncd_ssf_volume(im_n,im_m);
+ end
+ mm = workList(wx).wl(ii,1) ;
+ nn = workList(wx).wl(ii,2) ;
+ dxx{wx}(ii,:) = [min(d1,d2),mm,nn];
+ end
+end
+
+idxDel = [];
+
+tblVelo = getStats(dxx,pData,idxDel,tblVFMD);
+[rho,pCorr] = corr((tblVelo.velocity),(tblVelo.dmd));
+rho2 = rho^2
+toc(kernelCorrStart)
+4;
diff --git a/src/MATLAB/ncdVSvfmd/regnetCorr.m b/src/MATLAB/ncdVSvfmd/regnetCorr.m
new file mode 100644
index 0000000000000000000000000000000000000000..cbcc894a9766bacd4c4c9eb7c89df7f9639736c9
--- /dev/null
+++ b/src/MATLAB/ncdVSvfmd/regnetCorr.m
@@ -0,0 +1,16 @@
+function [rmse1,rmse2] = regnetCorr(tblVelo)
+% tblVelo = getStats(dxx,pData,[],tblVFMD);
+% X = [tblVelo.velocity,tblVelo.md];
+t2 = tblVelo(tblVelo.dmd~=0,:);
+
+X1 = [t2.velocity];
+X2 = [t2.velocity,t2.md];
+T = t2.dmd;
+
+mdl1 = fitrnet(X1,T,'Activations','sigmoid','LayerSizes',[1]);
+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"))
+
diff --git a/src/MATLAB/ncdVSvfmd/vfmd_meta.mat b/src/MATLAB/ncdVSvfmd/vfmd_meta.mat
new file mode 100644
index 0000000000000000000000000000000000000000..9cffe29aed7df086bee08ddf99448c753115f65d
Binary files /dev/null and b/src/MATLAB/ncdVSvfmd/vfmd_meta.mat differ