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

imPreProcess w/ empty radii

parent f73c8c37
Branches
No related tags found
No related merge requests found
function imf = getImf(pData,radii,scanType) function imf = getImf(pData,radii)
% radii is 1-element or 3-element cell array. [0,0,0] for unused channel % radii is 1-element or 3-element cell array. [0,0,0] for unused channel
% FLIF takes 1 or 3 UBYTE % FLIF takes 1 or 3 UBYTE
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 = {[0.5,8,8]}; % plate
%
% NGPU = HIP.Cuda.DeviceCount; % NGPU = HIP.Cuda.DeviceCount;
% p = ljsStartParallel(8*NGPU); % p = ljsStartParallel(8*NGPU);
...@@ -19,7 +20,7 @@ while startTarget <= size(pData,1) ...@@ -19,7 +20,7 @@ while startTarget <= size(pData,1)
if spmdIndex <= N_CUDA_PROCS && pp <= size(pData,1) if spmdIndex <= N_CUDA_PROCS && pp <= size(pData,1)
strDB = fullfile(pData.folder{pp},pData.filename{pp}); strDB = fullfile(pData.folder{pp},pData.filename{pp});
[im] = SSF.loadImage(strDB,1); [im] = SSF.loadImage(strDB,1);
% 3-channel plate,blob,0 % 3-channel e.g. plate,blob,0
for c = 1:length(radii) for c = 1:length(radii)
if all(0==radii{c}) if all(0==radii{c})
im1(:,:,:,c) = zeros(size(im)); im1(:,:,:,c) = zeros(size(im));
...@@ -28,6 +29,9 @@ while startTarget <= size(pData,1) ...@@ -28,6 +29,9 @@ while startTarget <= size(pData,1)
end end
end end
end end
% NGPU = HIP.Cuda.DeviceCount;
% p = ljsStartParallel(8*NGPU);
end end
startTarget = startTarget + N_CUDA_PROCS; startTarget = startTarget + N_CUDA_PROCS;
for i = 1:N_CUDA_PROCS for i = 1:N_CUDA_PROCS
...@@ -37,52 +41,52 @@ while startTarget <= size(pData,1) ...@@ -37,52 +41,52 @@ while startTarget <= size(pData,1)
end end
end end
% px = {}; % critical voxels px = {}; % critical voxels
% mx = []; % mean by stack mx = []; % mean by stack
% sx = []; % sigma by stack sx = []; % sigma by stack
% for i = 1 : length(imf) for i = 1 : length(imf)
% for c = 1:size(imf{i},4) for c = 1:size(imf{i},4)
% i1 = imf{i}(:,:,:,c); i1 = imf{i}(:,:,:,c);
% px{i,c} = i1(find(i1)); px{i,c} = i1(find(i1));
% mx(i,c) = mean(px{i,c}); mx(i,c) = mean(px{i,c});
% sx(i,c) = std(px{i,c}); sx(i,c) = std(px{i,c});
% end end
% end end
% clipRange = [2.5,97.5]; clipRange = [2.5,97.5];
% % clipRange = [0.5,99.5]; % clipRange = [0.5,99.5];
% for c = 1:size(imf{1},4) for c = 1:size(imf{1},4)
% p1 = vertcat(px{:,c}); p1 = vertcat(px{:,c});
% px_pos = p1(p1>0); px_pos = p1(p1>0);
% px_neg = abs(p1(p1<0)); px_neg = abs(p1(p1<0));
%
% clip_pos(c,:) = prctile(px_pos,clipRange); clip_pos(c,:) = prctile(px_pos,clipRange);
% clip_neg(c,:) = prctile(px_neg,clipRange); clip_neg(c,:) = prctile(px_neg,clipRange);
% end end
% 4; 4;
% for i = 1:length(imf) for i = 1:length(imf)
% for c = 1 : size(imf{i},4) for c = 1 : size(imf{i},4)
% if all(0==imf{i}(:,:,c)) if all(0==imf{i}(:,:,c))
% continue continue
% end end
% imp = imf{i}(:,:,:,c); imp = imf{i}(:,:,:,c);
% imp(imp<0)=0; imp(imp<0)=0;
% imn = imf{i}(:,:,:,c); imn = imf{i}(:,:,:,c);
% imn(imn>0)=0; imn(imn>0)=0;
% imn=abs(imn); imn=abs(imn);
% imq_pos = SSF.quantize8(imp,clip_pos(c,:)); imq_pos = SSF.quantize8(imp,clip_pos(c,:));
% imq_neg = SSF.quantize8(imn,clip_neg(c,:)); imq_neg = SSF.quantize8(imn,clip_neg(c,:));
% if all(0==imq_pos) if all(0==imq_pos)
% imf{i}(:,:,:,c) = imq_neg; imf{i}(:,:,:,c) = imq_neg;
% else else
% % [0,127] for imq_neg, [128,255] imq_pos % [0,127] for imq_neg, [128,255] imq_pos
% imq_pos = (0.5 .* imq_pos) + uint8(128 .* (imq_pos>0)); imq_pos = (0.5 .* imq_pos) + uint8(128 .* (imq_pos>0));
% imq_neg = 0.5 .* imq_neg ; imq_neg = 0.5 .* imq_neg ;
% imf{i}(:,:,:,c) = (imq_neg + imq_pos); imf{i}(:,:,:,c) = (imq_neg + imq_pos);
% end end
%
% end end
%
% imf{i} = uint8(imf{i}); imf{i} = uint8(imf{i});
% end end
% 4; 4;
function [imf,bwMask] = imPreProcess(im,radii) function [imf,bw] = imPreProcess(im,radii)
% threshold mask
bw = imbinarize(im);
bw = bwareaopen(bw,1e3);
if isempty(radii)
% radii = {} ? return thresholded input image
im(~bw) = 0;
imf = im;
return
end
im = uint8(im); % values are already uint8, just wrapped in a double im = uint8(im); % values are already uint8, just wrapped in a double
AX = 1/sqrt(3); AX = 1/sqrt(3);
......
...@@ -29,7 +29,10 @@ end ...@@ -29,7 +29,10 @@ end
tblMeta(tblMeta.bDisqualified,:) = []; tblMeta(tblMeta.bDisqualified,:) = [];
p = ljsStartParallel(96); p = ljsStartParallel(96);
imf = RSF.getImf(tblMeta,[[0.5,8,8] ; [5,5,5]],''); radii = {[0.5,8,8],[5,5,5],[0,0,0]}; % plate,blob
% radii = {[0.5,8,8]}; % plate
imf = RSF.getImf(tblMeta,radii);
d = []; d = [];
parfor i = 1:length(imf) parfor i = 1:length(imf)
dj = []; dj = [];
......
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
tic
ROOT = '/g/leverjs/Schuman_OCT/OCT/03-22-2024';
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');
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);
imf = RSF.getImf(tblMeta,[[0.5,8,8] ; [5,5,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,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;
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')
legend({'OS','OD'})
xlabel('NCD1')
ylabel('NCD2')
zlabel('NCD3')
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
\ No newline at end of file
load('onh_kernel_corr_pt5_5_15.mat') % load('onh_kernel_corr_pt5_5_15.mat')
% load('kcorr_VFMD_plate_blob.mat') % load('kcorr_VFMD_plate_blob.mat')
% load('kernel_corr_raw.mat') load('kernel_corr_raw.mat')
tblQualify = readtable('nyu onh matched.xlsx'); tblQualify = readtable('nyu onh matched.xlsx');
mx = regexp(tblQualify.file_name,'P(\d+)_(.+?)_(\d+-\d+-\d+)_(\d+-\d+-\d+)_(\w\w)_(.+?)_.*.img','tokens'); 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.scanID = cellfun(@(x) x{1}{6},mx,'UniformOutput',false);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment