goLR.m

tic
ROOT = '/g/leverjs/Schuman_OCT/OCT/qualifiedYes';
% ROOT = '/g/leverjs/Schuman_OCT/OCT/qualified';

PID = 'P10010';
% PID = 'P10125'
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}';
% database query
tblMeta = tblMeta(contains(tblMeta.scanType,'Optic'),:);
% tblMeta = tblMeta(contains(tblMeta.scanType,'Macula'),:);
% tblMeta = tblMeta(contains(tblMeta.eye,'OD'),:);

p = ljsStartParallel(96);

% radii = {[0.5,8,8],[5,5,5],[0,0,0]}; % plate,blob
% 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 = {[5,5,5],[8]}; % plate,blob, dark tube
radii

imf = RSF.getImf(tblMeta,radii);
% imf = RSF.quantizeImX(imf,[0,100]);
imf = RSF.quantizeImX(imf,[2.5,97.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,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')
zlabel('time')
legend({'OD','OS'})
xlabel('NCV1')
ylabel('NCV2')  
zlabel('time (months)')
title(['PID = ' PID ' :: normalized compression vectors (radii = '  jsonencode(radii) ')'])
% compute CSF
csf_mean = []; csf_std = [];
for k = 1:10
    [idx,Y] = Cluster.SpectralCluster(A,k);
    csf = CSF.csf_spatial(Y,idx);
    csf_mean(k) = mean(csf);
    csf_std(k) = std(csf);
end
% plot CSF
figure(2)
clf;errorbar(csf_mean,csf_std);xlim([0,10])
xlabel('number of clusters (K)')
ylabel('optimality deficiency')
title(['PID = ' PID ' :: CSF(radii = '  jsonencode(radii) ')'])

4;