diff --git a/readme.md b/readme.md
new file mode 100644
index 0000000000000000000000000000000000000000..848a065f75927f9b8061934f8ec3224387d00a50
--- /dev/null
+++ b/readme.md
@@ -0,0 +1,6 @@
+This is a repository for the manuscript and accompanying code currently under development.
+
+Once the manuscript is released, all code will be licensed open source (license TBD, likely GPL or BSD).
+
+questions? contact andrew.r.cohen 'at' drexel.edu
+
diff --git a/src/+RSF/checkMetaTiming.m b/src/+RSF/checkMetaTiming.m
new file mode 100644
index 0000000000000000000000000000000000000000..615c46b819d91b0db12e5527640e7671cd347abf
--- /dev/null
+++ b/src/+RSF/checkMetaTiming.m
@@ -0,0 +1,23 @@
+% function [pData] = getMetadata(ROOT,N,scanType,eye)
+scanType = 'Macula';
+eye = 'OD';
+tInterval = 0.5; % target time (years) for OCT image samples
+ROOT = '/g/leverjs/Schuman_OCT/09-21-2023';
+metaData = readtable(fullfile(ROOT,'09-21-2023 Data.xlsx'));
+[~,ic,ia] = unique(metaData.group);
+metaData.groupID = ia;
+
+[C,ia,ic] = unique(metaData.subjectID);
+N=[0,0,0];
+deltaT = [];
+for pp = 1:length(C)
+    pid = C(pp);
+    p1 = getPatientMetadata(metaData,pid,scanType,eye,tInterval);
+    if isempty(p1)
+        continue
+    end
+    dt = years(diff(unique(p1.dx)));
+    deltaT = [deltaT;dt];
+    N(p1.groupID(1)) = N(p1.groupID(1))+1; 
+end
+    
\ No newline at end of file
diff --git a/src/+RSF/chooseImageSamples.m b/src/+RSF/chooseImageSamples.m
new file mode 100644
index 0000000000000000000000000000000000000000..2ca9eb9bf771ccf54c291b8998963deb49fde299
--- /dev/null
+++ b/src/+RSF/chooseImageSamples.m
@@ -0,0 +1,26 @@
+function outData = chooseImageSamples(patientData,md1,tInterval)
+sampleDates = years(patientData.dx - min(patientData.dx));
+dt = unique(sampleDates);
+
+cxx = [];
+for i = 1:length(dt)
+    for j = i+1:length(dt)
+        for k = j+1:length(dt)
+            c1 = [dt(j)-dt(i), dt(k)-dt(j)];
+            c2 = vecnorm(c1-tInterval);
+            cxx = [cxx; i,j,k,c1,c2];
+        end
+    end
+end
+[~,idx] = min(cxx(:,end));
+idxT = cxx(idx,1:3);
+targetT = dt(idxT);
+outData = patientData(ismember(sampleDates,dt(idxT)),:);
+
+dt = years(diff(unique(outData.dx)));
+% dt is the interval between visits
+if any(dt>1) || any(dt < 1/12)
+    outData = [];
+end
+
+4;
diff --git a/src/+RSF/clipRSF.m b/src/+RSF/clipRSF.m
new file mode 100644
index 0000000000000000000000000000000000000000..49a0375e311cb604204c36b07c04aa97ea34c497
--- /dev/null
+++ b/src/+RSF/clipRSF.m
@@ -0,0 +1,15 @@
+function imf = clipRSF(imf)
+
+clipFactor = 20;
+clipDim = round(size(imf) ./ clipFactor);
+
+% clip long edge
+imf(1:clipDim(1),:,:)=0;
+imf(end-clipDim(1):end,:,:)=0;
+
+% clip two short edges
+imf(:,1:clipDim(2),:)=0;
+imf(:,end-clipDim(2):end,:)=0;
+
+imf(:,:,1:clipDim(3))=0;
+imf(:,:,end-clipDim(3):end)=0;
diff --git a/src/+RSF/cohortOCT.m b/src/+RSF/cohortOCT.m
new file mode 100644
index 0000000000000000000000000000000000000000..9154b3977b8797b13c2bcfe61f5772db9d2cd386
--- /dev/null
+++ b/src/+RSF/cohortOCT.m
@@ -0,0 +1,77 @@
+function tblResults = cohortOCT(ROOT,tblMeta,tblMD,nIterations,n_per_group,tInterval,scanType,radii)
+octStart = tic();
+
+tblResults = table();
+tblVelo = table();
+
+for nRep = 1:nIterations
+
+    pData = RSF.getMetadata(tblMeta,tblMD, n_per_group,scanType,'OD',tInterval);
+
+    imf = RSF.getImf(ROOT,pData,radii,scanType);
+    p = ljsStartParallel();
+    nPool = p.NumWorkers;
+    L = height(pData);
+    d = zeros(L,L);
+    cmdList = NCD.dParallelCommandList(d,nPool);
+    W = size(cmdList,2);
+    H = size(cmdList,1);
+    dxx = [];
+
+    parfor i = 1:H
+        mPrev = -1;
+        d1 = [];
+        for j = 1:W
+            m = cmdList(i,j,2);
+            n = cmdList(i,j,1);
+            if 0 == m || 0 == n
+                continue
+            end
+
+            im_m = imf{m};
+            im_n = imf{n};
+            d1 = SSF.ncd_ssf_volume(im_m,im_n);
+            if m==n
+                d2 = d1;
+            else
+                d2 = SSF.ncd_ssf_volume(im_n,im_m);
+            end
+            dxx(i,j) = min(d1,d2);
+        end
+    end
+
+    d = [];
+    for i = 1 : size(dxx,1)
+        for j = 1 : size(dxx,2)
+            c = cmdList(i,j,1);
+            r = cmdList(i,j,2);
+            if 0 == r || 0 == c
+                continue
+            end
+            d(r,c) = dxx(i,j);
+        end
+    end
+
+    for i = 1 : size(d,1)
+        for j = 1 : size(d,2)
+            if 0 == d(i,j)
+                d(i,j) = d(j,i);
+            end
+        end
+    end
+
+    tElapsedRep = toc(octStart)
+
+    [tv,pAnova] = RSF.progressionStats(pData,d); 
+    nt = table();
+    nt.pData = {pData};
+    nt.d = {d};
+    nt.pAnova = pAnova;
+    nt.tv = {tv};
+    tblResults = [tblResults;nt];
+    
+end
+
+
+
+
diff --git a/src/+RSF/corrOCT.m b/src/+RSF/corrOCT.m
new file mode 100644
index 0000000000000000000000000000000000000000..086da33452797fb8063462775eb5a8d4c2b88903
--- /dev/null
+++ b/src/+RSF/corrOCT.m
@@ -0,0 +1,76 @@
+function tblResults = corrOCT(pData,radii,scanType)
+octStart = tic();
+
+tblResults = table();
+tblVelo = table();
+
+imf = RSF.getImf(pData,radii,scanType);
+
+p = ljsStartParallel(64);
+nPool = p.NumWorkers;
+L = height(pData);
+d = zeros(L,L);
+cmdList = NCD.dParallelCommandList(d,nPool);
+W = size(cmdList,2);
+H = size(cmdList,1);
+dxx = [];
+
+parfor i = 1:H
+    mPrev = -1;
+    d1 = [];
+    for j = 1:W
+        m = cmdList(i,j,2);
+        n = cmdList(i,j,1);
+        if 0 == m || 0 == n
+            continue
+        end
+
+        im_m = imf{m};
+        im_n = imf{n};
+        d1 = SSF.ncd_ssf_volume(im_m,im_n);
+        if m==n
+            d2 = d1;
+        else
+            d2 = SSF.ncd_ssf_volume(im_n,im_m);
+        end
+        dxx(i,j) = min(d1,d2);
+    end
+end
+
+d = [];
+for i = 1 : size(dxx,1)
+    for j = 1 : size(dxx,2)
+        c = cmdList(i,j,1);
+        r = cmdList(i,j,2);
+        if 0 == r || 0 == c
+            continue
+        end
+        d(r,c) = dxx(i,j);
+    end
+end
+
+for i = 1 : size(d,1)
+    for j = 1 : size(d,2)
+        if 0 == d(i,j)
+            d(i,j) = d(j,i);
+        end
+    end
+end
+
+tElapsedRep = toc(octStart)
+
+% [tv,pAnova] = RSF.progressionStats(pData,d); 
+
+nt = table();
+nt.pData = {pData};
+nt.d = {d};
+nt.tblRes = {corrStatistics(pData,d)};
+% nt.pAnova = pAnova;
+% nt.tv = {tv};
+
+tblResults = [tblResults;nt];
+    
+
+
+
+
diff --git a/src/+RSF/draw1RSF.m b/src/+RSF/draw1RSF.m
new file mode 100644
index 0000000000000000000000000000000000000000..6feec6c550bd9a0d553afe2c466c463afa33cb1a
--- /dev/null
+++ b/src/+RSF/draw1RSF.m
@@ -0,0 +1,14 @@
+% pick the cohort for drawRSF()
+% load('iterOCT_10_20_30.mat')
+% load('iter_macula_10_20_30_N5.mat')
+radii = {10, 20, 30};
+ROOT = '/g/leverjs/Schuman_OCT/09-21-2023';
+scanType = 'Macula';
+% scanType = 'OpticNerve';
+
+idxTarget = 5; % pick 1 of tblResults cohorts. maybe from pAnova?
+d = tblResults.d{idxTarget};
+pAnova = tblResults.pAnova(idxTarget);
+pData = tblResults.pData{idxTarget};
+[~,Y] = Cluster.SpectralCluster(d,3);
+drawRSF(pData,pAnova,ROOT,radii,Y,scanType);
\ No newline at end of file
diff --git a/src/+RSF/drawOCT.m b/src/+RSF/drawOCT.m
new file mode 100644
index 0000000000000000000000000000000000000000..b4cfbb9a44b6efe3ba0cd7a25884ce43487da4d6
--- /dev/null
+++ b/src/+RSF/drawOCT.m
@@ -0,0 +1,35 @@
+
+A = d;
+K=3;
+[idx,Y] = Cluster.SpectralCluster(A,K);
+clf;hold off;
+    
+sym = {'ko','rx','m^','r*','gx','<k','b>','gs','bd'};
+
+
+[C,ia,ic] = unique(pData.group);
+
+hold on
+hx = [];
+for ff=1:size(Y,1)
+    idxSym = ic(ff);
+    h1 = plot3(Y(ff,1),Y(ff,2),Y(ff,3),sym{idxSym});
+    hx(idxSym) = h1;
+    
+%         id_str = strrep(tblOCT.PatientID{ff},'Healthy','H');
+%         id_str = strrep(id_str,'Glaucoma','G');
+%         id_str = strrep(id_str,'_M','');
+%         id_str = [id_str '_' tblOCT.AcquisitionDate{ff} '_' num2str(ff)];
+
+%         patient_id = tblOCT.subject_id(ff) - min(tblOCT.subject_id) +1;
+%         idx = find(tblOCT.subject_id == tblOCT.subject_id(ff));
+%         idxff = find(idx==ff);
+%         dx = datetime(tblOCT.date{ff});
+%         dx.Format = 'MMM-uuuu';
+%         id_str = [' ' num2str(patient_id) '_' char(dx) '_' tblOCT.eye{ff}];
+%         text(Y(ff,1),Y(ff,2),Y(ff,3),id_str,'Interpreter','none','color',color, 'FontSize',6)
+
+end
+set(gcf,'color','w')
+xlabel('u1');ylabel('u2'),zlabel('u3')
+legend(hx,C)
diff --git a/src/+RSF/drawPatients.m b/src/+RSF/drawPatients.m
new file mode 100644
index 0000000000000000000000000000000000000000..7a40b7616076daf27b9181ccfda5f31befe13b85
--- /dev/null
+++ b/src/+RSF/drawPatients.m
@@ -0,0 +1,40 @@
+[C,ia,ic] = unique(pData.group);
+[patient_C,patient_ia,patient_ic] = unique(pData.subjectID);
+patient_cmap = hsv(length(patient_C));
+A = d;
+K=3;
+[idx,Y] = Cluster.SpectralCluster(A,K);
+clf;hold off;
+
+sym = {'go','rx','b^','r*','gx','<k','b>','gs','bd'};
+hold on
+hx = [];
+
+for pp = 1:length(patient_C)
+    pidx = find(pData.subjectID == patient_C(pp));
+    idxSym = ic(pidx(1));
+    pxColor = sym{idxSym}(1);
+    pTbl = pData(pidx,:);
+    pY = Y(pidx,:);
+    pTime = pTbl.dx-min(pTbl.dx);
+    pY = [pY,years(pTime)];
+    visit = unique(pTime);
+    mxVisit = [];
+    for vv = 1:length(visit)
+        idxV = find(pTime == visit(vv));
+        yVisit = pY(idxV,:);
+        mxVisit(vv,:) = mean(yVisit,1);
+        for ss = 1:length(idxV)
+            plot3([yVisit(ss,1),mxVisit(vv,1)],[yVisit(ss,2),mxVisit(vv,2)],[yVisit(ss,end),mxVisit(vv,end)],':','color',pxColor,'LineWidth',0.1);
+        end
+        hx(idxSym) = plot3(mxVisit(vv,1),mxVisit(vv,2),mxVisit(vv,end),sym{idxSym});
+    end
+    for mm = 1:size(mxVisit,1)-1
+        plot3([mxVisit(mm,1),mxVisit(mm+1,1)],[mxVisit(mm,2),mxVisit(mm+1,2)],[mxVisit(mm,end),mxVisit(mm+1,end)],'-','color',pxColor)
+    end
+end
+
+set(gcf,'color','w')
+xlabel('k1');ylabel('k2'),zlabel('time (years)')
+legend(hx,C)
+axis square
\ No newline at end of file
diff --git a/src/+RSF/drawProgressionStats.m b/src/+RSF/drawProgressionStats.m
new file mode 100644
index 0000000000000000000000000000000000000000..9f7aeae4c2da4dabca5d548aa377c07b1240c537
--- /dev/null
+++ b/src/+RSF/drawProgressionStats.m
@@ -0,0 +1,18 @@
+
+
+% figure
+% boxplot(tblVelo.velocity,tblVelo.label,'notch',true)
+% set(gca,'YTickLabel',C)
+% ylabel('pattern velocity optic disc OCT on first 3 visits')
+
+
+% [p,tblAV,statsAV] = anova1(log(tblVelo.velocity),tblVelo.label);
+[p,tblAV,statsAV] = kruskalwallis(tblVelo.velocity,tblVelo.label);
+statsAV.gnames = {'-12';'-6:-12';'-6'};
+mc = multcompare(statsAV)
+
+% target = 'macula';
+target = 'optic disc';
+xlabel(['pattern velocity ' target ' OCT, 3 visits, ~6 month spacing'])
+title('')
+exportgraphics(gcf,['multcompare_' target '.png'])
diff --git a/src/+RSF/drawRSF.m b/src/+RSF/drawRSF.m
new file mode 100644
index 0000000000000000000000000000000000000000..ce9548b0e3bd4084a313ce9c2da9c73d0695e084
--- /dev/null
+++ b/src/+RSF/drawRSF.m
@@ -0,0 +1,126 @@
+% called by draw1RSF
+function drawRSF(pData,pAnova,ROOT,radius,Y,scanType)
+
+outdir = './rsfTiffs';
+if ~exist(outdir,'dir')
+    mkdir(outdir);
+end
+
+[C,ia,ic] = unique(pData.group);
+[patient_C,patient_ia,patient_ic] = unique(pData.subjectID,'rows','stable');
+
+imMontage = [];
+critPoints = {};
+top = 0;
+
+YY = {};
+rowMontage = {};
+NGPU = HIP.Cuda.DeviceCount;
+pPool = ljsStartParallel(4*NGPU);
+pix = {};
+imSeparate = {};
+parfor pp = 1:length(patient_C)
+    pidx = find(pData.subjectID == patient_C(pp));
+    pTbl = pData(pidx,:);
+    pTime = pTbl.dx-min(pTbl.dx);
+    visit = unique(pTime);
+    pTbl = pData(pidx,:);
+    pY = Y(pidx,:);
+    pY = [pY,years(pTime)];
+    imRow = [];
+    yVisit = [];
+    cp = {};
+    for vv = 1:length(visit)
+        idxV = find(pTime == visit(vv),1,'first');
+        strDB = fullfile(ROOT,scanType,pTbl.filename{idxV});
+        strDB = strrep(strDB,'.img','.LEVER');
+        im = SSF.loadImage(strDB,1);
+        [im1] = getImRSF(im,radius);
+        pix{pp}{vv} = im1(find(im1));
+        imSeparate{pp}{vv} = im1;
+        yVisit(vv,:) = pY(idxV,:);
+    end
+    YY{pp} = yVisit;
+    rowMontage{pp} = imRow;
+    critPoints{pp} = cp;
+end
+
+% cl = SSF.getClipQuantiles()
+pp = vertcat(pix{:});
+cl = SSF.getClipQuantiles(vertcat(pp{:}));
+imMontage = [];
+for pp = 1:length(imSeparate)
+    imRow = [];
+    for vv = 1:length(imSeparate{pp})
+        imSeparate{pp}{vv} = SSF.quantize8(imSeparate{pp}{vv},cl);
+        imRow = [imRow,imSeparate{pp}{vv}];
+    end
+    imMontage = [imMontage;imRow];
+end
+% imMontage = vertcat(rowMontage{:});
+figure(1);clf;imagesc(imMontage)
+hold on
+
+szim = [200,1024];
+xdim = szim(2);
+plot([xdim,xdim],[0,size(imMontage,1)],'-k','linewidth',2)
+plot(2.*[xdim,xdim],[0,size(imMontage,1)],'-k','linewidth',2)
+for i = 1:15
+    yy = (i-1)*szim(1);
+    if i==6
+        cc = 'r';
+    elseif i ==11
+        cc = 'm';
+    else
+        cc = 'k';
+    end
+
+    plot([0,size(imMontage,2)],[yy,yy],cc,'linewidth',2);
+end
+%     
+% for pp = 1:length(critPoints)
+%     
+%     ytop = (pp-1) * szim(1);
+%     for vv = 1:3
+%         xleft = (vv-1) * szim(2); 
+%         plot(xleft+critPoints{pp}{vv}(:,2),ytop+critPoints{pp}{vv}(:,1),'r.','markersize',8)
+%     end
+% end
+
+ylabels = {};
+ylocs = [szim(1)/2:szim(1):29*szim(1)/2];
+
+for pp = 1:length(patient_C)
+    gname = pData.group{patient_ia(pp)};
+    gname = strrep(gname,'Group ','')
+   ylabels{pp} = [num2str(patient_C(pp)) ' :: '  gname ];
+end
+set(gca,'YTick',ylocs)
+set(gca,'YTickLabel',ylabels);
+xinc = szim(2)/2:szim(2):5*szim(2)/2;
+set(gca,'XTick',xinc)
+set(gca,'XTickLabel',{'Visit 1','Visit 2','Visit 3'})
+
+for pp = 1:length(patient_C)
+    y = (pp-1)*szim(1) + szim(1)/2;
+    dy = diff(YY{pp});
+    v = vecnorm(dy(:,1:3),2,2) ./ dy(:,end);
+    text( szim(2),y,['V=' num2str(v(1),2) ',\Deltat=' num2str(dy(1,end),1)],...
+        'background','w','fontsize',6,'HorizontalAlignment','center','VerticalAlignment','middle')
+
+    text( 2 * szim(2),y,['V=' num2str(v(2),2) ',\Deltat=' num2str(dy(2,end),1)],...
+        'background','w','fontsize',6,'HorizontalAlignment','center','VerticalAlignment','middle')
+end
+
+cb = colorbar
+cb.Label.String = 'RSF';
+cb.Ticks = [1,64,128,192,255];
+cb.TickLabels = {'-1.0','-0.5','0.0','0.5','1.0'}
+tstr = [scanType ' :: pAnova = ' num2str(pAnova,2)];
+title(tstr);
+% fstr = ['p_' num2str(round(pAnova*100),2) '.png'];
+% exportgraphics(gcf,fullfile(outdir,fstr));
+exportgraphics(gcf,[scanType '_montage.pdf'])
+4;
+
+
diff --git a/src/+RSF/drawRSF_lever.m b/src/+RSF/drawRSF_lever.m
new file mode 100644
index 0000000000000000000000000000000000000000..98bd8bdd3485589ed3bfdcc380479e06926a0e29
--- /dev/null
+++ b/src/+RSF/drawRSF_lever.m
@@ -0,0 +1,35 @@
+% ROOT = '/g/leverjs/Schuman_OCT/09-21-2023/Macula';
+ROOT = '/g/leverjs/Schuman_OCT/09-21-2023/OpticNerve';
+outfolder = '/g/leverjs/Schuman_OCT/RSF_render_30';
+if ~exist(outfolder,'dir') 
+        mkdir(outfolder);
+end
+radii = 50; %{5,10,15,20};
+
+if iscell(radii)
+    rx = mat2str([radii{:}]);
+else
+    rx = mat2str(radii);
+end
+flist = dir(fullfile(ROOT,'**/*.LEVER'));
+
+for ff = 1:1 %length(flist)
+    strDB = fullfile(flist(ff).folder,flist(ff).name);
+    [im,CONSTANTS] = leversc.loadImage(strDB,1);
+    [imf,bwMask] = imPreProcess(im,radii);     
+    imc = mat2gray(LoG.separateLoG(imf));
+    imCombined = im;
+    imCombined(:,:,:,2:3) = imc;
+    imCombined(:,:,:,4) = bwMask;
+    CONSTANTS.imageData.NumberOfChannels = 4;
+    CONSTANTS.imageData.ChannelNames{1} = 'intensity image w/ NLM';
+    CONSTANTS.imageData.ChannelNames{2} = ['RSF LoG+ ' rx] ;
+    CONSTANTS.imageData.ChannelNames{3} = ['RSF LoG- ' rx] ;
+    CONSTANTS.imageData.ChannelNames{4} = ['bwMask'] ;
+
+    CONSTANTS.imageData.ChannelColors = [1,1,1;1,0,0;0,1,0];
+
+    MicroscopeData.WriterH5(imCombined,'path',outfolder,'imageData',CONSTANTS.imageData);
+end
+
+Import.leverImport('',outfolder)
\ No newline at end of file
diff --git a/src/+RSF/getDeltaTstats.m b/src/+RSF/getDeltaTstats.m
new file mode 100644
index 0000000000000000000000000000000000000000..6151f415d8fd644c8ee89e776723b4d40fb89ef9
--- /dev/null
+++ b/src/+RSF/getDeltaTstats.m
@@ -0,0 +1,19 @@
+% load('iter_OCT_10_20_30_N10_cp.mat')
+
+deltaT = [];
+labelT = [];
+for tt = 1:height(tblResults)
+    pData = tblResults.pData{tt};
+    [C,ia,ic] = unique(pData.subjectID);
+    for cc = 1:length(C)
+        patientC = find(pData.subjectID == C(cc));
+        patientData = pData(patientC,:);
+        dt = patientData.dx - min(patientData.dx);
+        dt = diff(unique(dt));
+        dt(dt == 0) = [];
+        deltaT = [deltaT;years(dt)];
+        labelT = [labelT;repmat(patientData.groupID(1),length(dt),1)];
+    end
+end
+[pAnova,tblAV,statsAV] = anova1(log(deltaT),labelT,'off');
+mc = multcompare(statsAV)
\ No newline at end of file
diff --git a/src/+RSF/getImRSF.m b/src/+RSF/getImRSF.m
new file mode 100644
index 0000000000000000000000000000000000000000..2606afe347f44849070ceb9322c24d5f01683fa7
--- /dev/null
+++ b/src/+RSF/getImRSF.m
@@ -0,0 +1,57 @@
+function [imrsf,cm] = getImRSF(im,sigma)
+
+% im = imresize3(im,[200,200,200]);
+% im is double, but on [0,255] (ish)
+% cast it to uint8 with no loss
+im = uint8(im);
+% 
+% if length(sigma) == 1
+%     sigma = [sigma, sigma, sigma];
+% end
+
+% if iscell(sigma)
+%     imp = single(0*im);
+%     imn = imp;
+%     for rr = 1:length(sigma)
+%         sx = sigma{rr};
+%         sx = [sx,sx,sx];
+%         cudaTarget = getCudaTarget;
+%         imf = HIP.LoG(im,sx,cudaTarget);
+%         [~,imp1,imn1] = LoG.separateLoG(imf);
+%         imp = max(imp,imp1);
+%         imn = max(imn,imn1);
+%     end
+%     imf = imp;
+%     imf(imn>imp) = -1.*imn(imn>imp);
+% else
+%     cudaTarget = getCudaTarget;
+%     imf = HIP.LoG(im,sigma,cudaTarget);
+% end
+
+imf = imPreProcess(im,sigma);
+
+
+imp = imf;
+imp(imp<0) = 0;
+imn = imf;
+imn(imn>0) = 0;
+imn = abs(imn);
+
+imp = imp .* imregionalmax(imp);
+imn = imn .* imregionalmax(imn);
+
+imp2 = max(imp,[],3);
+imn2 = max(imn,[],3);
+
+imrsf = imp2;
+imrsf(imn2>imp2) = -1 .* imn2(imn2>imp2);
+% imrsf = im2uint8(mat2gray(imrsf));
+imrsf = imrsf'; % put long axis horizontal for montage
+stepSize= (max(imrsf(:)) - min(imrsf(:))) ./ 255;
+% min + stepSize*nSteps0 = 0
+nSteps0 = ceil(-1 .* min(imrsf(:)) ./ (stepSize));
+cm = parula();
+cm(nSteps0+1,:) = [1,1,1];
+
+4;
+
diff --git a/src/+RSF/getImf.m b/src/+RSF/getImf.m
new file mode 100644
index 0000000000000000000000000000000000000000..d2d7887c4a7ae2d5a28571ba873f3b3a3e528e02
--- /dev/null
+++ b/src/+RSF/getImf.m
@@ -0,0 +1,40 @@
+function imf = getImf(pData,radii,scanType)
+
+NGPU = HIP.Cuda.DeviceCount;
+p = ljsStartParallel(8*NGPU);
+
+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;
+end
+
+px = px{:};
+px_pos = px(px>0);
+px_neg = abs(px(px<0));
+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)];
+
+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;
diff --git a/src/+RSF/getMetadata.m b/src/+RSF/getMetadata.m
new file mode 100644
index 0000000000000000000000000000000000000000..a199b7b751be4715ba199ec134924ab1746084a5
--- /dev/null
+++ b/src/+RSF/getMetadata.m
@@ -0,0 +1,30 @@
+function [pData] = getMetadata(metaData,tblMD,N,scanType,eye,tInterval)
+
+px = unique(metaData.subjectID);
+pid = randperm(length(px));
+
+pData = table();
+
+nFound = 0;
+nidx = 0;
+while(nFound < N && nidx < length(px) )
+    nidx = nidx + 1;
+    id1 = px(pid(nidx));
+    mdidx = find(tblMD.subjectID == id1 & strcmp(tblMD.Eye,eye));
+    if length(mdidx) < 2
+        continue
+    end
+    md1 = tblMD(mdidx,:);
+    md1.dx = datetime(md1.visit_date_date,'InputFormat','MM-dd-yyyy');
+    [C,ia,ic] = unique(md1.dx);
+    md1 = md1(ia,:);
+    p1 = RSF.getPatientMetadata(metaData,md1,id1,scanType,eye,tInterval);
+    if isempty(p1)
+        continue
+    end
+    p1.md = repmat(NaN,height(p1),1);
+    nFound = nFound + 1;
+    pData = [pData ; p1];
+end
+
+4;
diff --git a/src/+RSF/getPatientMetadata.m b/src/+RSF/getPatientMetadata.m
new file mode 100644
index 0000000000000000000000000000000000000000..6f4130d61866ecda0db2a08aece6caeba0100e7d
--- /dev/null
+++ b/src/+RSF/getPatientMetadata.m
@@ -0,0 +1,11 @@
+function pOut = getPatientMetadata(metaData,md1,pid,scanType,eye,tInterval)
+pOut = [];
+p1 = metaData(metaData.subjectID == pid,:); % current patient candidate
+
+p1 = sortrows(p1,'dx');
+d1 = unique(p1.dx);
+if length(d1) < 3
+    return
+end
+
+pOut = RSF.chooseImageSamples(p1,md1,tInterval);
diff --git a/src/+RSF/imPreProcess.m b/src/+RSF/imPreProcess.m
new file mode 100644
index 0000000000000000000000000000000000000000..0fbc9cdf46ea29f8164f7f8d1bba92f7d25a1b0d
--- /dev/null
+++ b/src/+RSF/imPreProcess.m
@@ -0,0 +1,43 @@
+function [imf,bwMask] = imPreProcess(im,radii)
+
+% AX = 1/sqrt(3) .* [15.38,1,1] ;
+% AX = 1/sqrt(3) .* [5,1,1] ;
+AX = 1/sqrt(3);
+% dx = 15.38;
+dx = 20;
+AX = 1/sqrt(3);
+% fun fact :: mat2gray here causes r2 decline
+im = uint8(im); % values are already uint8, just wrapped in a double
+
+imp = double(0*im);
+imn = imp;
+for rr = 1:length(radii)
+    sx = radii(rr);
+    if 1 == length(sx)
+        sx = [0.45,sx,sx];
+    end
+    sx = AX .* sx;
+    cudaTarget = getCudaTarget;
+    imf = HIP.LoG(im,sx,cudaTarget);
+    [~,imp1,imn1] = LoG.separateLoG(imf);
+    imp = max(imp,imp1);
+    imn = max(imn,imn1);
+end
+
+% threshold mask
+bw = imbinarize(im);
+bw = bwareaopen(bw,1e3);
+bw = imdilate(bw,ones(ceil(AX .* max(radii))));
+% ??? threshold before or after regionalmax ???
+% challenge problem -- connect this choice to K(X|M) 
+imn = imn .* imregionalmax(imn);
+imn = imn .* bw;
+imp = imp .* imregionalmax(imp);
+imp = imp .* bw;    
+
+imf = imp;
+imf(imn>=imp) = -1 * imn(imn>=imp);
+
+4;
+
+
diff --git a/src/+RSF/optimizeDeltaTstats.m b/src/+RSF/optimizeDeltaTstats.m
new file mode 100644
index 0000000000000000000000000000000000000000..9d67fb88b4e64624c9884c423ab3c010182eb2f4
--- /dev/null
+++ b/src/+RSF/optimizeDeltaTstats.m
@@ -0,0 +1,28 @@
+    
+ROOT = '/g/leverjs/Schuman_OCT/09-21-2023';
+metaData = readtable(fullfile(ROOT,'09-21-2023 Data.xlsx'));
+[px] = unique(metaData.subjectID);
+[G,ia,ic] = unique(metaData.group);
+metaData.groupID = ic;
+scanType = 'Macula'
+
+pData = table();
+pkw = [];
+for tInterval = 0.2:.1:2
+    deltaT = [];
+    labelT = [];
+    nskip = 0;
+    for pp = 1 : length(px)
+        p1 = getPatientMetadata(metaData,px(pp),scanType,'OD',tInterval);
+        if isempty(p1)
+            nskip = nskip + 1;
+            continue
+        end
+        dt = p1.dx - min(p1.dx);
+        dt = diff(unique(dt));
+        dt(dt == 0) = [];
+        deltaT = [deltaT;years(dt)];
+        labelT = [labelT;repmat(p1.groupID(1),length(dt),1)];
+    end
+    pkw = [pkw;tInterval,kruskalwallis(deltaT,labelT,'off'),nskip]
+end
diff --git a/src/+RSF/pickMinVelocityStack.m b/src/+RSF/pickMinVelocityStack.m
new file mode 100644
index 0000000000000000000000000000000000000000..faa820c5b18f4c9ade74156e340807bdb41b57b3
--- /dev/null
+++ b/src/+RSF/pickMinVelocityStack.m
@@ -0,0 +1,35 @@
+function [velo,md] = pickMinVelocityStack(pTbl,pY)
+
+pTime = pTbl.dx-min(pTbl.dx);
+pY = [pY,years(pTime)];
+[visit,iVisit,cVisit] = unique(pTime);
+assert(3==length(visit),'ASSERT FAILED # visits ~= 3. pickMinVelocityStack.m');
+md = pTbl.md(iVisit);
+md = md(2:end); % one for each velocity
+
+mxVisit = {};
+for vv = 1:length(visit)
+    mxVisit{vv} = [];
+    idxV = find(pTime == visit(vv));
+    for ii = 1:length(idxV)
+        yVisit = [pY(idxV(ii),:),vv,ii];
+        mxVisit{vv} = [mxVisit{vv};yVisit];
+    end
+end
+vxx = [];
+for i = 1:size(mxVisit{1},1)
+    for j = 1:size(mxVisit{2},1)
+        for k = 1:size(mxVisit{3},1)
+            vvv = [mxVisit{1}(i,1:end-2);mxVisit{2}(j,1:end-2);mxVisit{3}(k,1:end-2)];
+            dy = diff(vvv);
+            velo = vecnorm(dy(:,1:3),2,2) ./ dy(:,end);
+            vxx = [vxx;velo',i,j,k];
+        end
+    end
+end
+[~,idxMin] = min(vxx(:,1)+vxx(:,2));
+velo = vxx(idxMin,1:2)';
+
+
+
+
diff --git a/src/+RSF/pickRSFradius.m b/src/+RSF/pickRSFradius.m
new file mode 100644
index 0000000000000000000000000000000000000000..3d5a96adf83a26721607f42d037669c4aad0c975
--- /dev/null
+++ b/src/+RSF/pickRSFradius.m
@@ -0,0 +1,31 @@
+ROOT = '/g/leverjs/Schuman_OCT/09-21-2023';
+nIterations = 10;
+N_PER_GROUP = 5;
+% scanType = 'Macula';
+scanType = 'Optic Disc';
+tInterval = 0.5;
+
+pData = getMetadata(ROOT, N_PER_GROUP,scanType,'OD',tInterval);
+
+rx = [];
+
+cudaTarget = getCudaTarget;
+rx = repmat({[]},10,1);
+
+for pp = 1:height(pData)
+    if contains(scanType,'Optic')
+        strDB = fullfile(ROOT,'OpticNerve',pData.filename{pp});
+    else
+        strDB = fullfile(ROOT,scanType,pData.filename{pp});
+    end
+    strDB = strrep(strDB,'.img','.LEVER');
+
+    [im] = SSF.loadImage(strDB,1);
+    
+    for r = 1:10
+        sx = repmat(r,1,3);
+        imf = HIP.LoG(im,sx,cudaTarget);
+        rx{r} = [rx{r};pp,min(imf(:)),max(imf(:))];
+    end
+    cellfun(@(x) [x(end,1),mean(x(:,2:3),1)],rx,'UniformOutput',false)
+end
\ No newline at end of file
diff --git a/src/+RSF/platemetric.m b/src/+RSF/platemetric.m
new file mode 100644
index 0000000000000000000000000000000000000000..c371566c0c3c724974c5e2fd0c39cb6e3129730c
--- /dev/null
+++ b/src/+RSF/platemetric.m
@@ -0,0 +1,45 @@
+function imout = platemetric(imn,bwn)
+
+imn = mat2gray(imn);
+% mimic fibermetric, with thickness blur
+% thickness = 5;
+% sigma=thickness/6;
+% Ig = HIP.Gaussian(imn,[15.38,1,1],1,0);
+Ig = imn;
+
+% gSpace = round(size(imn)./min(size(imn)));
+gSpace = [15.38,1,1];
+% gSpace = [1,1,1];
+[gx,gy,gz] = gradient(Ig,gSpace(1),gSpace(2),gSpace(3));
+[gxx,gxy,gxz] = gradient(gx,gSpace(1),gSpace(2),gSpace(3));
+[gyx,gyy,gyz] = gradient(gy,gSpace(1),gSpace(2),gSpace(3));
+[gzx,gzy,gzz] = gradient(gz,gSpace(1),gSpace(2),gSpace(3));
+
+idx = find(bwn);
+imout = 0*imn;
+
+for i = 1:length(idx) %length(r)
+    idx1 = idx(i);    
+    Hyx = [ gxx(idx1), gxy(idx1), gxz(idx1);gyx(idx1),gyy(idx1),gyz(idx1);...
+        gzx(idx1),gzy(idx1),gzz(idx1)];
+    [evec,eigval] = eig(Hyx);
+    eigval = diag(eigval);
+    if all(0==eigval),continue,end
+    if any(isnan(eigval)),continue,end
+    [~,idxev] = sort(abs(eigval));
+    ev = eigval(idxev);
+    ev = ev ./ norm(ev);    
+    if ev(end)>0
+        continue
+    end
+    ev = abs(ev);
+    rb =  sqrt(ev(1)*ev(2)) / ev(3) ; % plate small
+    ra = ev(1) / ev(2); % plate small, tube large    
+
+    s = imn(idx1);
+%     v = (exp(-ra^2 / 0.5) ) * (exp(-rb^2 / 0.5)) * s(1 - exp(-s^2 / 2));
+    v =  (1-exp(-ra^2 / 0.5) ) * (exp(-rb^2 / 0.5)) * (1 - exp(-s^2 / 2));
+    imout(idx1) = v;
+end
+% imout = imout .* imbinarize(imout);
+4;
\ No newline at end of file
diff --git a/src/+RSF/progressionStats.m b/src/+RSF/progressionStats.m
new file mode 100644
index 0000000000000000000000000000000000000000..a5171f2ba0bbce55b569bf8d5647969ffbcd6805
--- /dev/null
+++ b/src/+RSF/progressionStats.m
@@ -0,0 +1,23 @@
+function [tblVelo,pAnova] = progressionStats(pData,d)
+% [C,ia,ic] = unique(pData.group);
+[patient_C,patient_ia,patient_ic] = unique(pData.subjectID);
+A = d;
+K=3;
+[idx,Y] = Cluster.SpectralCluster(A,K);
+
+tblVelo = table();
+for pp = 1:length(patient_C)
+    pidx = find(pData.subjectID == patient_C(pp));
+    pTbl = pData(pidx,:);
+    pY = Y(pidx,:);
+    [velo,md] = RSF.pickMinVelocityStack(pTbl,pY);
+    nt = table();
+    nt.velocity = (velo);
+    nt.md = (md);
+    tblVelo = [tblVelo;nt];
+end
+
+pAnova = 0;
+% [pAnova] = anova1(log(tblVelo.velocity),tblVelo.label,'off')
+% mc = multcompare(statsAV)
+
diff --git a/src/+RSF/thresholdRetinalOCT.m b/src/+RSF/thresholdRetinalOCT.m
new file mode 100644
index 0000000000000000000000000000000000000000..c058b903a582e62dd40a3c494af0a1727c3d2547
--- /dev/null
+++ b/src/+RSF/thresholdRetinalOCT.m
@@ -0,0 +1,9 @@
+function bw = thresholdRetinalOCT(im)
+
+bw = imbinarize(im);
+cc = bwconncomp(bw); 
+stats = regionprops3(cc,"Volume"); 
+idx = find([stats.Volume] >= 0.25 * max(stats.Volume) ); 
+bw = ismember(labelmatrix(cc),idx); 
+
+4;
\ No newline at end of file
diff --git a/src/readme_src.md b/src/readme_src.md
new file mode 100644
index 0000000000000000000000000000000000000000..b3b884a2ae197c243a4c3b5ca08aad95a2a4c3ff
--- /dev/null
+++ b/src/readme_src.md
@@ -0,0 +1,3 @@
+Sample applications and metadata manipulation TBD.
+
+Main repository for now on bioimage/leverApps/...contact AC for access.