From 53490aa3f698df28f5fbcf963ecfc678c2a3b845 Mon Sep 17 00:00:00 2001 From: Robert <52195853+robdahn@users.noreply.github.com> Date: Thu, 11 Jan 2024 16:48:03 +0100 Subject: [PATCH] Changed: Updated PBTsimple (see issue #19). Changed paths: M CHANGES.txt M cat_main_reportfig.m M cat_vol_pbtsimple.m --- CHANGES.txt | 8 +++ cat_main_reportfig.m | 6 +- cat_vol_pbtsimple.m | 163 +++++++++++++++++++++++++++++++++++++------ 3 files changed, 153 insertions(+), 24 deletions(-) diff --git a/CHANGES.txt b/CHANGES.txt index 651e5a20..ba6d83c2 100644 --- a/CHANGES.txt +++ b/CHANGES.txt @@ -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: diff --git a/cat_main_reportfig.m b/cat_main_reportfig.m index 672ce80d..f0ab5246 100644 --- a/cat_main_reportfig.m +++ b/cat_main_reportfig.m @@ -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 % -------------------------------------------------------------- @@ -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') diff --git a/cat_vol_pbtsimple.m b/cat_vol_pbtsimple.m index 815cb2a6..66eee1f9 100644 --- a/cat_vol_pbtsimple.m +++ b/cat_vol_pbtsimple.m @@ -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); + 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 @@ -129,11 +214,15 @@ % 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; @@ -141,7 +230,34 @@ % -------------------------------------------------------------------- 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; % 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 @@ -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 @@ -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; @@ -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. % ...