Skip to content
Snippets Groups Projects
Commit c6d61f5b authored by ac's avatar ac
Browse files

rsf.darktubes rnfl et al forte

parent 1b472bc6
Branches
No related tags found
No related merge requests found
function imfm = darkTubes(im,CONSTANTS,thickness)
voxelSize = CONSTANTS.imageData.PixelPhysicalSize;
voxelSize([1,2]) = voxelSize([2,1]); % swap (x,y,z) to matlab column major (y,x,z)
vx = max(voxelSize)./voxelSize;
% vx = [1,1,1];
sigma = 1/sqrt(3) * thickness * vx;
cudaTarget = getCudaTarget;
imf = HIP.LoG(im,sigma,cudaTarget);
imn = imf;
imn(imn>0) = 0;
imn = abs(imn);
bw = imbinarize(im);
bw = bwareaopen(bw,1e3);
bw = imdilate(bw,ones(4*thickness));
imn(~bw) = 0;
imfm = fibermetric(imn);
% imfm = imfm .* imregionalmax(imfm);
4;
......@@ -19,16 +19,22 @@ while startTarget <= size(pData,1)
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);
% RSF images and metadata via lever file
[im,CONSTANTS] = leversc.loadImage(strDB,1,1);
% 3-channel e.g. plate,blob,0
for c = 1:length(radii)
if all(0==radii{c})
im1(:,:,:,c) = zeros(size(im));
else
if length(radii{c}) == 1
% hardcode for darktubes...careful...ACK
im1(:,:,:,c) = RSF.darkTubes(im,CONSTANTS,radii{c});
else
im1(:,:,:,c) = RSF.imPreProcess(im,radii{c});
end
end
end
end
% NGPU = HIP.Cuda.DeviceCount;
% p = ljsStartParallel(8*NGPU);
......
% clipLimits = [minVal,maxVal];
function [i1q,i2q] = quantizePair(i1,i2)
if isa(i1,'uint8') && isa(i2,'uint8')
return
end
% if isa(i1,'uint8') && isa(i2,'uint8')
% return
% end
pix = [i1(find(i1));i2(find(i2))];
......
tic
ROOT = '/g/leverjs/Schuman_OCT/OCT/combined';
ROOT = '/g/leverjs/Schuman_OCT/OCT/qualified';
PID = 'P10010';
flist = dir(fullfile(ROOT,['/' PID '*.LEVER']));
......@@ -13,25 +13,12 @@ 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);
radii = {[0.5,8,8],[5,5,5],[0,0,0]}; % plate,blob
% 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)
......
function iRNFL = getRNFL(p1,tblRNFL)
% find p1 in tblRNFL and return avg thickness
[~,f1,~] = fileparts(p1.filename);
[~,tblFname,~] = fileparts(tblRNFL.file_name);
idx1 = find(strcmp(f1,tblFname),1,'first');
iRNFL = tblRNFL{idx1,"AVERAGETHICKNESS"};
4;
function tblVelo = getStats(dxx, pData, tblMD)
if ~exist('idxDel','var')
idxDel = [];
function tblVelo = getStats(dxx, pData, tblMD, tblRNFL)
if ~exist('tblRNFL','var')
tblRNFL = [];
end
tblVelo = table();
idxUnused = 1:height(pData);
for i = 1:length(dxx)
for j = 1:size(dxx{i},1)
jRNFL = getRNFL(pData(j,:),tblRNFL);
m = dxx{i}(j,2);
n = dxx{i}(j,3);
% m and n are indices into pData
% check first if either is disqualified
if pData.disqualified(m) || pData.disqualified(n)
continue
end
if pData.dx(m) > pData.dx(n)
4;
......@@ -30,26 +32,37 @@ for i = 1:length(dxx)
if dmd>0
continue
end
mRNFL = getRNFL(pData(m,:),tblRNFL);
nRNFL = getRNFL(pData(n,:),tblRNFL);
deltaRNFL = nRNFL - mRNFL;
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.ncd = (dxx{i}(j,1));
nt.md = min(pData.md(n),pData.md(m));
nt.dmd = dmd;
nt.mn = [m,n];
nt.ddx = pData.dx(n) - pData.dx(m);
nt.mRNFL = mRNFL;
nt.nRNFL = nRNFL;
nt.deltaRNFL = deltaRNFL;
% find sex
idx = find(tblMD.subjectID == pData.subjectID(m),1,'first');
nt.sex = tblMD.gender_full(idx);
if abs(nt.dmd)>4
continue
end
% on the unqualified data vs. qualified data
% if years(nt.ddx)>1
% continue
% end
% if nt.velocity > 0.9 && nt.velocity < 0.975
% if nt.ncd > 0.9 && nt.ncd< 0.975
% continue
% end
if m ~= n
......
......@@ -10,7 +10,7 @@ fx = {'nyu VF qualified data great than -6 data export.xlsx',...
tblVFMD = table();
for ff = 1 : length(fx)
tx = readtable(fullfile("importMetadata/dataExcel/",fx{ff}));
tx = readtable(fullfile("../importMetadata/dataExcel/",fx{ff}));
% filter 1
idxVFMD = contains(lower(tx.qualification_status),'yes');
tx = tx(idxVFMD,:);
......
datetime
kernelCorrStart = tic();
eye = 'OD';
ROOT = '/g/leverjs/Schuman_OCT/OCT/combined';
ROOT = '/g/leverjs/Schuman_OCT/OCT/qualified';
% 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
% tblVFMD = getVFMD(eye);
tblMeta = getMetaOCT(ROOT,eye,scanType);
load('../qualify/onh_rnfl_thickness.mat');
radii = {[0.5,8,8],[5,5,5],[4]}; % plate,blob
% radii = {[0.5,8,8]}; % plate
[rho2,tblVelo, pData, dxx] = kernelCorr(radii,tblVFMD,tblMeta,tblRNFL);
[rmse1,~] = regnetCorr(tblVelo);median(rmse1)
function [rho2,tblVelo, pData, dxx] = kernelCorr(radii,tblVFMD,tblMeta,scanType,eye)
function [rho2,tblVelo, pData, dxx] = kernelCorr(radii,tblVFMD,tblMeta,tblRNFL)
kernelCorrStart = tic();
......@@ -40,7 +40,7 @@ end
p=ljsStartParallel(96);
imf = RSF.getImf(pData,radii,scanType);
imf = RSF.getImf(pData,radii);
[C,ia,ic] = unique(pData.subjectID);
......@@ -84,8 +84,8 @@ parfor wx = 1:length(workList)
end
end
tblVelo = getStats(dxx,pData,tblVFMD);
[rho,pCorr] = corr((tblVelo.velocity),(tblVelo.dmd));
tblVelo = getStats(dxx,pData,tblVFMD,tblRNFL);
[rho,pCorr] = corr((tblVelo.ncd),(tblVelo.dmd));
rho2 = rho^2
toc(kernelCorrStart)
4;
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];
% X = [tblVelo.ncd,tblVelo.md];
t2 = tblVelo(tblVelo.dmd~=0,:);
X1 = [t2.velocity];
X2 = [t2.velocity,t2.md];
X1 = [t2.ncd];
X2 = [t2.ncd,t2.md];
T = t2.dmd;
% T = t2.deltaRNFL;
% leave one out x-validation
err = [];
mdlX = {};
xTarget = X1;
xTarget = X2;
for i = 1:length(xTarget)
trainX = xTarget;
trainX(i,:) = [];
......@@ -27,8 +27,8 @@ for i = 1:length(xTarget)
mdlX{i} = mdl;
pred = mdl.predict(xTarget(i,:));
err(i) = abs(T(i) - pred);
[i,mean(err)]
[i,median(err)]
end
mean(err)
% mean(err)
4;
% getRNFL
tx_onh = readtable('nyu onh matched.xlsx');
save('onh_rnfl_thickness.mat','tx_onh');
ROOT = '/g/leverjs/Schuman_OCT/OCT/combined';
target = '/g/leverjs/Schuman_OCT/OCT/qualified';
if ~exist(target,'dir')
mkdir(target);
end
% flist = dir(fullfile(ROOT,'**/*.LEVER'));
% qONH = readtable('nyu onh matched.xlsx');
% qMacula = readtable('nyu mac matched.xlsx');
% [~,fONH,~] = fileparts(qONH.file_name);
% [~,fMacula,~] = fileparts(qMacula.file_name);
qlist = [];
for i = 1 : length(flist)
[~,fname,~] = fileparts(flist(i).name);
idx = find(strcmp(fname,fMacula));
if ~isempty(idx)
tidx = qMacula(idx,:);
else
idx = find(strcmp(fname,fONH));
tidx = qONH(idx,:);
end
if isempty(idx)
continue
end
qs = lower(tidx.qualification_status);
if any(strcmp(qs,'no'))
continue
end
qlist = [qlist,flist(i)];
src = fullfile(flist(i).folder,fname);
cmd = ['cp "' src '".* ' target];
system(cmd);
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment