Skip to content

Commit

Permalink
Changed: Updated PBTsimple (see issue #19).
Browse files Browse the repository at this point in the history 
Changed paths:
 M CHANGES.txt
 M cat_main_reportfig.m
 M cat_vol_pbtsimple.m
  • Loading branch information
@robdahn
robdahn committed Jan 11, 2024
1 parent 7ad2f2f commit 53490aa
Show file tree
Hide file tree
Showing 3 changed files with 153 additions and 24 deletions.
8 changes: 8 additions & 0 deletions CHANGES.txt
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,12 @@
------------------------------------------------------------------------
r2518 | dahnke | 2024-01-11 16:47:58

Changed paths:
M CHANGES.txt
M cat_main_reportfig.m
M cat_vol_pbtsimple.m

Changed: Updated PBTsimple (see issue #19).------------------------------------------------------------------------
r2517 | gaser | 2024-01-09 16:15:15

Changed paths:
Expand Down
6 changes: 3 additions & 3 deletions cat_main_reportfig.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1358,7 +1358,7 @@ function cat_main_reportfig(Ym,Yp0,Yl1,Psurf,job,qa,res,str)
if job.extopts.print>1
if exist('Psurf','var') && ~isempty(Psurf)
if 1 %~strcmpi(spm_check_version,'octave') && opengl('info')
boxwidth = 0.2;
boxwidth = 0.1;
if job.extopts.report.type <= 1
%% classic top view
% --------------------------------------------------------------
Expand Down Expand Up @@ -1483,10 +1483,10 @@ function cat_main_reportfig(Ym,Yp0,Yl1,Psurf,job,qa,res,str)
end
maxdiff = 4 * ceil(std(cdata(:))*8)/8;
srange = [-maxdiff maxdiff];
boxwidth = diff(srange)/40; % 0.1;
boxwidth = diff(srange)/40 / 2; % 0.05;
else
srange = [0 6];
boxwidth = diff(srange)/30; % 0.2;
boxwidth = diff(srange)/30 / 2; % 0.1;
end
%% hrange = srange(1) + boxwidth/2:boxwidth:srange(2);
if strcmpi(renderer,'opengl')
Expand Down
163 changes: 142 additions & 21 deletions cat_vol_pbtsimple.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,14 +32,99 @@

if ~exist('opt','var'), opt = struct(); end

def.levels = 4; % 0 .. 4
def.refine = 1; % refined WM distance based on CSF distance (myelination correction)
def.gyrusrecon = 1; % additional PBT gyri reconstruction
def.correctoffeset = 2;
def.extendedrange = 1;
def.WMTC = 1; % Partial voxel-based topology correction of small
% WM-defects based on the volume output of cat_vol_genus0
def.levels = 4; % Number of dual distance estimate (requires refinement) with 0
% for the classic approach with 1.5 and 2.5 as CSF and WM boudnary.
% Good values are between 2 and 4, whease higher values (>6)
% run into problems for interpolation artifacts of close to the
% full tissue class.
def.refine = 1; % Refined WM distance based on CSF distance (myelination correction)
def.gyrusrecon = 1; % additional PBT gyri reconstruction (TRUE/FALSE)
% this reduce thickness overestimations but maybe underestimates
% the outer surface position in gyral regions (running to much inside)
def.correctoffeset = 2; % not really required if no refinement is used
def.extendedrange = 1; % may causes closed gyri
def.sharpening = 1; % sharpening the Ypp map to support more better resampling to lower resolution for the 0.5 layer - worse
opt = cat_io_checkinopt(opt,def);


% == preparations ==
% should maybe be better part of create CS

if 0
% INITIAL SHARPENING?
% RD202401: remove this part when the rest is finished
% not so powerfull as you use already the multple distance estimates
% should also not work in a label map!

% sharpen WM
Yp0x = min(3,max(2,Yp0));
Yp0x = min(3,max(2,Yp0x + .5*(Yp0x - cat_vol_smooth3X(Yp0x,1)) + .5*(Yp0x - cat_vol_smooth3X(Yp0x,2)) ));
Yp0(Yp0(:)>2) = Yp0x(Yp0(:)>2);

% sharpen CSF
Yp0x = min(2,max(1,Yp0));
Yp0x = min(2,max(1,Yp0x + .5*(Yp0x - cat_vol_smooth3X(Yp0x,1)) + .5*(Yp0x - cat_vol_smooth3X(Yp0x,2)) ));
Yp0(Yp0(:)<2 & Yp0(:)>1) = Yp0x(Yp0(:)<2 & Yp0(:)>1);
end


if 1 % opt.WMTC
% RD20231224: for WM topology smoothing
% the idea is to apply the voxel-based correction and _close_ small WM holes
% (<10 voxel) for the 2.75 and 2.25 boundaries for values above 2.

% indicate isolated holes and replace by median of the neighbors
Yp0(Yp0<0.35 & ~cat_vol_morph(Yp0<1,'l')) = 1; % close major wholes in the WM
Ymsk = Yp0==-1 & cat_vol_morph(Yp0>0.9,'d',1); % filter small wholes close to the WM
Yp0 = cat_vol_median3(single(Yp0)+1,Ymsk,~Ymsk)-1;

% indicate isolated objects and replace by median of the neighbors
Yp0(Yp0>0.65 & cat_vol_morph(Yp0==-1,'l')) = -1;
Ymsk = Yp0>0.95 & cat_vol_morph(Yp0<-0.9,'d',1);
Yp0 = cat_vol_median3(single(Yp0)+1,Ymsk,~Ymsk)-1;

% RD20210401: but there are some new background dots
Ymsk = Yp0==-1 & cat_vol_morph(Yp0>-1,'lc'); % close major wholes in the WM
Yp0 = cat_vol_median3(single(Yp0)+1,Ymsk,~Ymsk)-1;
end



if opt.WMTC
% topology correction for WM based on the volume output of cat_vol_genus0

% Closing of the WM object to avoid smaller holes.
% Opening is here not useful as the WM !
tclevels = [2.75, 2.25];
for ii = 1:1
for tci = 1:numel(tclevels)
evalc(sprintf('Yppc = cat_vol_genus0(Yp0, %0.2f,0);',tclevels(tci)));
%if tclevels(tci) > 2.5
Yholes = Yppc==1 & Yp0<tclevels(tci);
Yholes = Yholes & ~cat_vol_morph(Yholes,'l',[1e4,8]); % allow only small corrections
Yp0(Yholes) = tclevels(tci) + 0.01; % close
clear Yholes;
%else
Ybridges = Yppc==0 & Yp0>=tclevels(tci);
Ybridges = Ybridges & ~cat_vol_morph(Ybridges,'l',[1e4,8]); % allow only small corrections
Yp0(Ybridges) = tclevels(tci) - 0.01; % open
clear Ybridges;
%end
end
end


% object correction in the background (opening of GM-WM objects)
tclevels = [1.75, 1.25, 1.01];
for tci = 1:numel(tclevels)
Yp0(Yp0>tclevels(tci) & ~cat_vol_morph(Yp0>tclevels(tci), ...
'ldo',2.5 - tclevels(tci))) = tclevels(tci) - 0.01;
end
end


% (1) Distance estimation:
% --------------------------------------------------------------------
if opt.levels == 0
Expand Down Expand Up @@ -129,19 +214,50 @@
% If gyri were reconstructed too than also the WMD would have to be
% corrected to avoid underestimation of the position map with surfaces
% running to close to the WM.
overrange = .0; % range extention of the Ypp was resulting in worse
% T1 and position values and more topology defects
YM = Yp0 > 1.5-0.45*opt.extendedrange & Yp0 < 2.5+0.45*opt.extendedrange & Ygmt>eps;
Ycdc = Ycd; Ycdc(YM) = min(Ycd(YM), Ygmt(YM) - Ywd(YM));
Ypp = zeros(size(Yp0),'single'); Ypp(Yp0 >= 2.5+0.45*opt.extendedrange) = 1;
Ypp(YM) = Ycdc(YM) ./ (Ygmt(YM) + eps);
Ypp(Ypp>1) = 0;
Yppl = -2 + 2*smooth3(Ypp>0); Yppl(Ypp>0) = 0;
Ypph = 2*smooth3(Ypp>=1); Ypph(Ypp<1) = 0;
Ypp = min(1 + overrange,max(0 - overrange, Ypp + Yppl + Ypph ));
clear Ycdc YM;


% (5) Voxel-size resolution correction:
% --------------------------------------------------------------------
Ygmt = Ygmt * mean(vx_vol);


if 0
% FINAL TOPOLOGY CORRECTION
% final voxel-based topology correction of some levels
% however, the topology correction in the surface create should be more powerfull
tclevels = 1:-0.25:0;
for tci = 1:numel(tclevels)
evalc(sprintf('Yppc = cat_vol_genus0(Ypp, %0.2f,0);',tclevels(tci)));
if tclevels(tci) > 0.5
Ypp(Yppc==1 & Ypp<tclevels(tci)) = tclevels(tci) + 0.01; % close
else
Ypp(Yppc==0 & Ypp>=tclevels(tci)) = tclevels(tci) - 0.01; % open
end
end
end


% sharpening - important to compensate the downsampling
if opt.sharpening
Ypp = Ypp + 1*(Ypp - cat_vol_smooth3X(Ypp,1/mean(vx_vol))) + 2*(Ypp - cat_vol_smooth3X(Ypp,2/mean(vx_vol))) + 4*(Ypp - cat_vol_smooth3X(Ypp,4/mean(vx_vol)));

% final smoothing to prepare surface reconstruction that also correct for WM topology issues
spm_smooth(Ypp,Ypp,.5/mean(vx_vol));
Ypp = min(1 + overrange,max(0 - overrange,Ypp));

end
end
% ======================================================================
function [Ycd, Ywd] = cat_vol_cwdist(Yp0,opt)

% CSF and WM distance
Expand All @@ -151,14 +267,15 @@
% multi-level distance estimation
hss = opt.levels; % number of opt.levels (as pairs)
for si = 1:hss
offset = max(0,min(.5,si * .5/hss)) / 2;
offset = max(0,min(.5, si * (.5)/hss)) / 2;

Ycdl = cat_vbdist(single(Yp0 < ( 1.5 - offset)), Yp0 < 2.5 + 0.45*opt.extendedrange ); Ycdl(Ycdl > 1000) = 0;
Ycdh = cat_vbdist(single(Yp0 < ( 1.5 + offset)), Yp0 < 2.5 + 0.45*opt.extendedrange ); Ycdh(Ycdh > 1000) = 0;

Ycdl = cat_vbdist(single(Yp0 < ( 1.5 - offset)), Yp0 < 2.5 + 0.45*opt.extendedrange );
Ycdh = cat_vbdist(single(Yp0 < ( 1.5 + offset)), Yp0 < 2.5 + 0.45*opt.extendedrange );

if opt.extendedrange
Ycdlc = max(0,cat_vbdist(single(Yp0 < ( 2.5 - offset)), Yp0 < 2.95 ) - (opt.levels<8));
Ycdhc = max(0,cat_vbdist(single(Yp0 < ( 2.5 + offset)), Yp0 < 2.95 ) - (opt.levels<8));
% additiona correction for levels>8
Ycdlc = max(0,cat_vbdist(single(Yp0 < ( 2.5 - offset)), Ycdl>0 ) - .5 ); Ycdlc(Ycdlc > 1000) = 0;
Ycdhc = max(0,cat_vbdist(single(Yp0 < ( 2.5 + offset)), Ycdh>0 ) - .5 ); Ycdhc(Ycdhc > 1000) = 0;
Ycdl = Ycdl - Ycdlc;
Ycdh = Ycdh - Ycdhc;
end
Expand All @@ -172,21 +289,24 @@
Ycdl(Ycdl>0) = max(eps, Ycdl(Ycdl>0) - (.5 - offsetc ));
Ycdh(Ycdh>0) = max(eps, Ycdh(Ycdh>0) + (.5 - offsetc ));
end
Ycd = Ycd + .5/hss .* Ycdl + .5/hss .* Ycdh;
Ycw = max(0.1,.5/hss - opt.refine * Ycd/10 );

% mixing
Ycd = Ycd + .5/hss .* Ycdl + .5/hss .* Ycdh;
Ycw = max(0.1,.5/hss - opt.refine * Ycd/5 );
%%
clear Ycdl Ycdh;


% WM distances
Ywdl = cat_vbdist(single(Yp0 > ( 2.5 - offset)), Yp0 > 1.5 - 0.45*opt.extendedrange );
Ywdh = cat_vbdist(single(Yp0 > ( 2.5 + offset)), Yp0 > 1.5 - 0.45*opt.extendedrange);
Ywdl = cat_vbdist(single(Yp0 > ( 2.5 - offset)), Yp0 > 1.5 - 0.45*opt.extendedrange ); Ywdl(Ywdl > 1000) = 0;
Ywdh = cat_vbdist(single(Yp0 > ( 2.5 + offset)), Yp0 > 1.5 - 0.45*opt.extendedrange ); Ywdh(Ywdh > 1000) = 0;

if opt.extendedrange
Ywdlc = max(0,cat_vbdist(single(Yp0 > ( 1.5 - offset)), Yp0 > 1.05 ) - (opt.levels<8));
Ywdhc = max(0,cat_vbdist(single(Yp0 > ( 1.5 + offset)), Yp0 > 1.05 ) - (opt.levels<8));
Ywdlc = max(0,cat_vbdist(single(Yp0 > ( 1.5 - offset)), Yp0 > 1.05 ) - .5); Ywdlc(Ywdlc > 1000) = 0;
Ywdhc = max(0,cat_vbdist(single(Yp0 > ( 1.5 + offset)), Yp0 > 1.05 ) - .5); Ywdhc(Ywdhc > 1000) = 0;
Ywdl = Ywdl - Ywdlc;
Ywdh = Ywdh - Ywdhc;
end

if opt.correctoffeset
if opt.correctoffeset==2
offsetc = median(Ywdl(Ywdl>0) - Ywdh(Ywdl>0))/2;
Expand All @@ -196,13 +316,14 @@
Ywdl(Ywdl>0) = max(eps, Ywdl(Ywdl>0) + (.5 - offsetc ));
Ywdh(Ywdh>0) = max(eps, Ywdh(Ywdh>0) - (.5 - offsetc ));
end

% mixing
Ywd = Ywd + Ycw .* Ywdl + (.5/hss*2 - Ycw) .* Ywdh;
end

Ycd(Ycd>100) = 0;
Ywd(Ywd>100) = 0;

end
% ======================================================================
function Yd = cat_vol_eudist(Yb,Ymsk,levels,correctoffeset)
%cat_vol_eudist. Euclidean distance estimation to mixed boundaries.
% ...
Expand Down

0 comments on commit 53490aa

Please sign in to comment.