goLR.m

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

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}';

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 = {[3]} % dark tubes
imf = RSF.getImf(tblMeta,radii);
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;
    idxOD = find(strcmp(tblMeta.eye,'OD'));
    idxOS = find(strcmp(tblMeta.eye,'OS'));
    clf;hold on
    if size(Y,2) > 2
        plot3(Y(idxOD,1),Y(idxOD,2),Y(idxOD,3),'r*')
        plot3(Y(idxOS,1),Y(idxOS,2),Y(idxOS,3),'og')
        zlabel('NCD3')
    else
        plot(Y(idxOD,1),Y(idxOD,2),'r*')
        plot(Y(idxOS,1),Y(idxOS,2),'og')
        
    end
    legend({'OD','OS'})
    xlabel('NCD1')
    ylabel('NCD2')
    
% 
% 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

idxOD = find(strcmp(tblMeta.eye,'OS'));
aOD = A(idxOD,idxOD);

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
clf;errorbar(csf_mean,csf_std);xlim([0,10])