diff --git a/cat_conf_opts.m b/cat_conf_opts.m index 7b567442..896f4c45 100644 --- a/cat_conf_opts.m +++ b/cat_conf_opts.m @@ -39,14 +39,6 @@ }; tpm.num = [1 1]; tpm.def = @(val)cat_get_defaults('opts.tpm', val{:}); -if expert>1 - tpm.help = [tpm.help - { - 'In case of multiple TPMs, the original image is affine registrated to each TPM and a kmeans algorithm is used to fit 3 classes. The best two images with the lowest standard deviations are combined to a new less biased TPM. The TPMs must include the same classes in MNI space!' - '' - }]; - tpm.num = [1 inf]; -end % ngaus: diff --git a/cat_main.m b/cat_main.m index e9f12eb9..dce21d2f 100644 --- a/cat_main.m +++ b/cat_main.m @@ -265,7 +265,7 @@ if job.extopts.LASstr>1 || job.extopts.inv_weighting [Ymi,Ymt,Ycls] = cat_main_LASs(Ysrc,Ycls,Ym,Yb,Yy,Tth,res,vx_vol,extoptsLAS2); % use Yclsi after cat_vol_partvol else - [Ymi,Ymt,Ycls] = cat_main_LAS2(Ysrc,Ycls,Ym,Yb,Yy,T3th,res,vx_vol,job.extopts,Tth); + [Ymi,Ymt,Ycls] = cat_main_LAS(Ysrc,Ycls,Ym,Yb,Yy,T3th,res,vx_vol,job.extopts,Tth); end clear Ymt; % use original global scaled ... stime2 = clock; % not really correct but better than before diff --git a/cat_main1639.m b/cat_main1639.m index f595dbe5..c1a6986e 100644 --- a/cat_main1639.m +++ b/cat_main1639.m @@ -298,7 +298,7 @@ if job.extopts.LASstr>1 [Ymi,Ym,Ycls] = cat_main_LASs(Ysrc,Ycls,Ym,Yb,Yy,Tth,res,vx_vol,extoptsLAS2); % use Yclsi after cat_vol_partvol else - [Ymi,Ym,Ycls] = cat_main_LAS2(Ysrc,Ycls,Ym,Yb,Yy,T3th,res,vx_vol,job.extopts,Tth); + [Ymi,Ym,Ycls] = cat_main_LAS(Ysrc,Ycls,Ym,Yb,Yy,T3th,res,vx_vol,job.extopts,Tth); end stime2 = clock; % not really correct but better than before diff --git a/cat_main_LAS.m b/cat_main_LAS.m index f30e4798..3d813243 100644 --- a/cat_main_LAS.m +++ b/cat_main_LAS.m @@ -1,120 +1,155 @@ -function [Yml,Ymg,Ycls,Ycls2,T3th] = cat_main_LAS(Ysrc,Ycls,Ym,Yb0,Yy,T3th,res,vx_vol,extopts,Tth) +function [Yml,Ymg,Ycls,Ycls2,T3th] = ... + cat_main_LAS(Ysrc,Ycls,Ym,Yb0,Yy,T3th,res,vx_vol,extopts,Tth) % This is an exclusive subfunction of cat_main. % ______________________________________________________________________ % % Local Adaptive Segmentation (LAS): % % This version of the local adaptive intensity correction includes a -% bias correction that is based on a maximum filter for the WM and a mean -% filter of GM to stabilize the correction in region with less WM. -% The extension is mostly based on the assumption that the tissue next to -% the CSF (and high divergence sulci) has to be WM (maximum, high -% divergence) or GM. For each tissue a refined logical map is generated -% and used to estimate the local intensity threshold. -% It is important to avoid high intensity blood vessels in the process, -% because they will push down local WM and GM intensity - due to the CSF -% near possition of blood vessels that mostly push down GM. -% Based on these values an intensity transformation is used. Compared to -% the global correction this has to be done for each voxel. To save time -% only a rough linear transformation is used. -% -% Finally, a second NLM-filter is used and a refinement of WM structures -% by a divergence map +% bias correction that is based on a maximum filter for the WM, a mean +% filter of GM and a mean filter of the head tissue (if available).m +% For each tissue a refined logical map is generated to filter the local +% tissue intensity and approximate the local tissue intensity in the whole +% volume. Based on these values an intensity transformation is used. +% Compared to the global correction this has to be done for each voxel. +% To save time only a rough linear transformation is used. Finally, a +% second NLM-filter is used. +% It is important to avoid high intensity blood vessels in the process +% because they will push down local WM and GM intensity. +% % ______________________________________________________________________ % -% [Yml,Ycls,Ycls2,T3th] = ... -% cat_main_LAS(Ysrc,Ycls,Ym,Yb0,Yy,T3th,res,vx_vol,PA,template) +% [Yml,Ymg,Ycls,Ycls2,T3th] = ... +% cat_main_LAS(Ysrc,Ycls,Ym,Yb0,Yy,T3th,res,vx_vol,extopts,Tth) % -% Yml .. local intensity correct image -% Ycls .. corrected SPM tissue class map -% Ycls2 .. ? -% T3th .. tissue thresholds of CSF, GM, and WM in Ysrc +% Yml .. local intensity normalized image +% Ymg .. global intensity normalized image +% Ycls .. corrected SPM tissue class map +% Ycls2 .. ? +% T3th .. tissue thresholds of CSF, GM, and WM in Ysrc +% +% Ysrc .. (bias corrected) T1 image +% Ycls .. SPM tissue class map +% Ym .. intensity corrected T1 image (BG=0,CSF=1/3,GM=2/3,WM=1) +% Yb0 .. brain mask +% Yy .. deformation map +% T3th .. intensity thresholds from global intensity normalization +% res .. SPM segmentation structure +% vx_vol .. voxel dimensions +% extopts .. +% Tth .. +% ______________________________________________________________________ % -% Ysrc .. (bias corrected) T1 image -% Ym .. intensity corrected T1 image (BG=0,CSF=1/3,GM=2/3,WM=1) -% Ycls .. SPM tissue class map -% Yb0 .. brain mask -% Yy .. deformation map -% res .. SPM segmentation structure -% vx_vol .. voxel dimensions -% PA .. CAT atlas map -% template .. ? +% Christian Gaser, Robert Dahnke +% Structural Brain Mapping Group (https://neuro-jena.github.io) +% Departments of Neurology and Psychiatry +% Jena University Hospital +% ______________________________________________________________________ +% $Id$ + + + % ______________________________________________________________________ % % internal maps: % % Yg .. gradient map - edges between tissues -% Ydiv .. divergence map - sulci, gyris pattern, and blood vessels +% Ydiv .. divergence map - sulci/gyri pattern, and blood vessels % Yp0 .. label map - tissue classes (BG=0,CSF=1,GM=2,WM=3) % % Ysw .. save WM tissue map % Ybv .. blood vessel map % Ycp .. CSF / background area for distances estimation -% Ycd .. CSF / background distance +% Ycd .. CSF / background distance map +% Ywd .. WM distance map +% % Ycm .. CSF % Ygm .. GM % Ywm .. WM % Yvt .. WM next to the ventricle map % ______________________________________________________________________ % -% Christian Gaser, Robert Dahnke -% Structural Brain Mapping Group (https://neuro-jena.github.io) -% Departments of Neurology and Psychiatry -% Jena University Hospital +% Development / TODO: +% * replace old calls of cat_vol_morph by distance based operations with +% resolution parameter % ______________________________________________________________________ -% $Id$ - def.uhrlim = 0.7; - extopts = cat_io_checkinopt(extopts,def); - extopts.uhrlim = 0.2; % the reduction to 1 mm was not realy good for the Ycls map + + % The reduction to 1 mm is not really good for the Ycls map. If this + % is required to support faster and memory saving preprocessing also + % for ultra-high resolution data further test are necessary! + %def.uhrlim = 0.7; + %extopts = cat_io_checkinopt(extopts,def); + extopts.uhrlim = 0.2; % no reduction for >0.4 mm - % set this variable to 1 for simpler debuging without reduceBrain + % set this variable to 1 for simpler debugging without reduceBrain % function (that normally save half of processing time) - verb = extopts.verb-1; - vxv = 1/ mean(vx_vol); - dsize = size(Ysrc); - NS = @(Ys,s) Ys==s | Ys==s+1; % function to ignore brain hemisphere coding - LASstr = max(eps,min(1,extopts.LASstr)); % LAS strenght (for GM/WM threshold)3 - manual correction based on R1109 (2017/02) - LAB = extopts.LAB; % atlas labels - LABl1 = 1; % use atlas map - cleanupstr = min(1,max(0,extopts.gcutstr)); % required to avoid critical regions + verb = extopts.verb - 1; + vxv = 1 / mean(vx_vol); % normalization of voxel size (mostly for old calls of cat_vol_morph) + dsize = size(Ysrc); + NS = @(Ys,s) Ys==s | Ys==s+1; % function to ignore brain hemisphere coding + LASstr = max(eps,min(1,extopts.LASstr)); % LAS strength (for GM/WM threshold) - manual correction based on R1109 (2017/02) + LAB = extopts.LAB; % atlas labels + LABl1 = 1; % use atlas map + cleanupstr = min(1,max(0,extopts.gcutstr)); % required to avoid critical regions (only used in case of atlas maps) cleanupdist = min(3,max(1,1 + 2*cleanupstr)); + % set debug = 1 and do not clear temporary variables if there is a breakpoint in this file - dbs = dbstatus; debug = 0; - for dbsi=1:numel(dbs), if strcmp(dbs(dbsi).name,'cat_main_LAS'); debug = 1; break; end; end + dbs = dbstatus; debug = 0; + for dbsi = 1:numel(dbs), if strcmp(dbs(dbsi).name,'cat_main_LAS'); debug = 1; break; end; end + + %% --------------------------------------------------------------------- % First, we have to optimize the segments using further information that % SPM do not use, such as the gradient, divergence and distance maps. -% The gradient map (average of the first derivate of the T1 map) is an -% edge map and independent of the image intensity. It helps to avoid PVE -% regions and meninges. -% The divergence (second derivate of the T1 map) help to identfiy sulcal -% and gyral pattern and therefore to find WM and CSF regions for furhter -% corrections and to avoid meninges and blood vessels. -% Furhtermore, special assumption can be used. -% The first one is the maximum property of the WM in T1 data that allows -% using of a maxim filter for the GM/WM region. +% (1) The gradient map (average of the first derivate of the T1 map) is +% an edge map and independent of the image intensity. It helps to +% avoid PVE regions and meninges. +% (2) The divergence (second derivate of the T1 map) help to identify +% sulcal and gyral pattern and therefore to find WM and CSF regions +% for further corrections and to avoid meninges and blood vessels. +% (3) Distance maps can help to describe the cortex by its expected +% thickness or to separate structures close to the skull or deep +% in the brain. +% +% Furthermore, special assumption can be used. +% The first is that the maximum property of the WM in T1 data that +% allows using of a maxim filter for the GM/WM region. % The second is the relative stable estimation of CSF/BG that allows to -% estimat a distance map. Because, most regions have a thin layer of -% GM around the WM we can avoid overestimation of the WM by the other +% estimate a distance map. Because, most regions have a thin layer of +% GM around the WM, we can avoid overestimation of the WM by the other % maps (especially the divergence). % --------------------------------------------------------------------- - fprintf('\n'); - stime = cat_io_cmd(' Prepare maps','g5','',verb); dispc=1; + fprintf('\n'); + stime = cat_io_cmd(' Prepare maps','g5','',verb); + + + % general resolution limitation + % ------------------------------------------------------------------- + % Using a lot of maps make this function memory intensive that is + % critical for high resolution data. Furthermore, it can be expected + % that the full resolution is not required here. However, reducing + % the tissue class map can lead to other problems. Hence, the uhrlim + % parameter is maybe inactive (see above). + % ------------------------------------------------------------------- if any( vx_vol < extopts.uhrlim/2 ) - Yclso2 = Ycls; Ysrco2 = Ysrc; + % store old data that is needed in full resolution + Yclso2 = Ycls; + Ysrco2 = Ysrc; + + % reduce all input maps (Ysrc, Ycls, Yy, Ym, Yb0) [Ysrc,resT0] = cat_vol_resize( Ysrc , 'reduceV' , vx_vol , extopts.uhrlim , 64 ); - for ci=1:6 + for ci = 1:numel(Ycls) Ycls{ci} = cat_vol_resize( Ycls{ci}, 'reduceV' , vx_vol , extopts.uhrlim , 64 ); end Yy2 = ones([size(Ysrc),4],'single'); - for ci=1:3 + for ci = 1:numel(Yy) Yy2(:,:,:,ci) = cat_vol_resize( Yy(:,:,:,ci) , 'reduceV' , vx_vol , extopts.uhrlim , 64 ); end Yy = Yy2; clear Yy2; @@ -127,25 +162,35 @@ end - - % brain segmentation can be restricted to the brain to save time + % brain bounding box to save processing time Yclso = Ycls; - [Ym,BB] = cat_vol_resize( Ym , 'reduceBrain' , vx_vol , round(10/mean(vx_vol)) , Yb0 ); - Yb = cat_vol_resize( Yb0 , 'reduceBrain' , vx_vol , round(10/mean(vx_vol)) , Yb0 ); clear Yb0; - for i=1:6, Ycls{i} = cat_vol_resize(Ycls{i} , 'reduceBrain' , vx_vol , BB.BB); end + [Ym,BB] = cat_vol_resize( Ym , 'reduceBrain' , vx_vol , round(10/mean(vx_vol)) , Yb0 ); + Yb = cat_vol_resize( Yb0 , 'reduceBrain' , vx_vol , round(10/mean(vx_vol)) , Yb0 ); clear Yb0; + try + for i = 1:numel(Ycls), Ycls{i} = cat_vol_resize(Ycls{i} , 'reduceBrain' , vx_vol , BB.BB); end + catch + error('LAS failed probably due to bad tissue contrast.\n'); + end - % helping maps (Yg = mean gradient = edge) and divergence + % helping maps (Yg = mean gradient = edge) and divergence (CSF or WM skeleton) Yg = cat_vol_grad( Ym , vx_vol ); Ydiv = cat_vol_div( max(0.33,Ym) , vx_vol ); - Yp0 = single( Ycls{1} )/255*2 + single( Ycls{2} )/255*3 + single( Ycls{3} )/255; - Yb = smooth3( Yb | cat_vol_morph(Yb,'d',2*vxv) & Ym<0.8 & Yg<0.3 & Ym>0 )>0.5; % increase brain mask, for missing GM + Yp0 = single( Ycls{1} )/255*2 + single( Ycls{2} )/255*3 + single( Ycls{3} )/255; % create label map + Yb = smooth3( Yb | cat_vol_morph(Yb,'d',2*vxv) & Ym<0.8 & Yg<0.3 & Ym>0 )>0.5; % increase brain mask, for missing GM + + - %% adding of atlas information (for subcortical structures) + %% adding of atlas and template information (for subcortical structures) + % ------------------------------------------------------------------- +%%% 20180801 - Why not use the partitioning before LAS to benefit by the + % better ventricle, WMH and cerebellum maps? + % - Because the partitioning benefit in subcortical regions + % from LAS. >> call it twice? % ------------------------------------------------------------------- if LABl1 - stime = cat_io_cmd(' Prepare partitions','g5','',verb,stime); dispc=dispc+1; + stime = cat_io_cmd(' Prepare partitions','g5','',verb,stime); % map atlas to RAW space for i=1:5 @@ -153,168 +198,307 @@ Vl1A = spm_vol(extopts.cat12atlas{1}); break catch - % read error in parallel processing - pause(1) + % avoid read error in parallel processing + pause(rand(1)) end end - Yl1 = cat_vol_ctype(round(spm_sample_vol(Vl1A,double(Yy(:,:,:,1)),double(Yy(:,:,:,2)),double(Yy(:,:,:,3)),0))); - Yl1 = reshape(Yl1,dsize); + Yl1 = cat_vol_ctype( cat_vol_sample(res.tpm(1),Vl1A,Yy,0) ); + % load WM of the Dartel/Shooting Template for WMHs - use uint8 to save memory + % Ywtpm .. WM template map Vtemplate = spm_vol(extopts.templates{end}); - Ywtpm = cat_vol_ctype(spm_sample_vol(Vtemplate(2),double(Yy(:,:,:,1)),double(Yy(:,:,:,2)),double(Yy(:,:,:,3)),0)*255,'uint8'); - if ~debug, clear Yy; end - Ywtpm = single(reshape(Ywtpm,dsize)); spm_smooth(Ywtpm,Ywtpm,2*vxv); Ywtpm = cat_vol_ctype(Ywtpm); + Ywtpm = cat_vol_ctype( cat_vol_sample(res.tpm(1),Vtemplate(2),Yy,1) * 255,'uint8'); + if debug==0, clear Yy; end + Ywtpm = single(reshape(Ywtpm,dsize)); + spm_smooth(Ywtpm,Ywtpm,2*vxv); + Ywtpm = cat_vol_ctype(Ywtpm); + + + % apply boundary box for brain mask Yl1 = cat_vol_resize(Yl1 ,'reduceBrain',vx_vol,round(4/mean(vx_vol)),BB.BB); Ywtpm = cat_vol_resize(Ywtpm,'reduceBrain',vx_vol,round(4/mean(vx_vol)),BB.BB); + % The correction of image resolution is not required because of the adaptation of the Yy? + - LASmod = min(2,max(1,mean((Ym( NS(Yl1,LAB.BG) & Yg<0.1 & Ydiv>-0.05 & Ycls{1}>4)) - 2/3) * 8)); % do not reduce LASstr + % do not reduce LASstr + LASmod = min(2,max(1,mean((Ym( NS(Yl1,LAB.BG) & Yg<0.1 & Ydiv>-0.05 & Ycls{1}>4)) - 2/3) * 8)); + else + LASmod = 1; %#ok end - %% adaptation of the LASstr depending on average basal values - LASstr = min(1,max(0.01,LASstr * LASmod)); % adaptation by local BG variation - LASfs = 1 / max(0.01,LASstr); % smoothing filter strength - LASi = min(8,round(LASfs)); % smoothing iteration (limited) + % adaptation of the LASstr depending on average basal values + LASstr = min(1,max(0.01,LASstr * LASmod)); % adaptation by local BG variation + LASfs = 1 / max(0.01,LASstr); % smoothing filter strength + LASi = min(8,round(LASfs)); % smoothing iteration (limited) + + %% helping segments - stime = cat_io_cmd(sprintf(' Prepare segments (LASmod = %0.2f)',LASmod),'g5','',verb,stime); dispc=dispc+1; - % Ybb = don't trust SPM too much by using Yp0 because it may miss some areas! Shood be better now with MRF. + % ------------------------------------------------------------------- + % We will now try to identify the CSF and WM regions those finally + % describe the GM area. + % * brain mask (Ybb), blood vessel maps (Ybv), + % * distance maps from CSF (Ycd), WM (Ywd), and brain mask (Ybd) + % * tissue maps for CSF (Ycp >> Ycm), GM (Ygm), and WM (Ysw >> Ywm) + % * ventricle map (Yvt) + % ------------------------------------------------------------------- + stime = cat_io_cmd(sprintf(' Prepare segments (LASmod = %0.2f)',LASmod),'g5','',verb,stime); + + % Ybb .. brain mask + % Don't trust SPM too much by using Yp0 because it may miss some GM areas! + % 20180801 - The brain mask should be correct now. Ybb = cat_vol_morph((Yb & Ym>1.5/3 & Ydiv<0.05) | Yp0>1.5,'lo',vxv); - % Ysw = save WM and blood vessels mpas - % Ybv = possible blood vessels - Ysw = cat_vol_morph(Ycls{2}>128 & (min(1,Ym)-Ydiv)<1.5,'lc',vxv*2) & (Ym-Ydiv)>5/6; % 1.2 - Ybv = ((min(1,Ym) - Ydiv + Yg)>2.0 | (Ycls{5}>16 & Ym<0.6 & Ycls{1}<192)) & ... - ~cat_vol_morph(Ysw,'d',1) & Ym>0.2; - % Ycp = for CSF/BG distance initialization - Ycp = (Ycls{3}>240 & Ydiv>0 & Yp0<1.1 & Ym<0.5) | ... % typcial CSF - (Ycls{5}>8 & Ycls{2}<32 & Ym<0.6 & Ydiv>0) | ... % venes - ((Ym-Ydiv/4<0.4) & Ycls{3}>4 & Ycls{3}>16) | ... % sulcal CSF - ...(single(Ycls{6})+single(Ycls{5})+single(Ycls{4}))>192 | ... % save non-csf - (single(Ycls{6})+single(Ycls{4}))>192 | ... % save non-csf .. class 5 with error for negative t1 values - ~cat_vol_morph(Ybb,'lc',5) | ... % add background - Ym<0.3; % but do not trust the brain mask! - Ywd = cat_vbdist(single(Yp0>2.5),Yp0>0.5,vx_vol); % WM distance for skelelton - %Ysk = cat_vol_div(min(5,Ywd),2); %clear Ywd; % divergence skeleton - %Ysk = (Ym + min(0,Ysk))<0.2; % binary divergence skeleton - %Ycp = Ycp | Ysk; %clear Ysk; % - Ycp(smooth3(Ycp)>0.4)=1; % remove some meninges - Ycd = cat_vbdist(single(Ycp),~Ycp,vx_vol); % real CSF distance - Ycd((Ym-Ydiv<2/3 | Ydiv>0.1) & Ycls{3}>4 & Ycls{3}>1) = ... correction for sulci ... maybe a second distance estimation??= - min(Ycd((Ym-Ydiv<2/3 | Ydiv>0.1) & Ycls{3}>4 & Ycls{3}>1),1.5); - % we need to remove strong edge regions, because here is no GM layer between CSF and WM ??? - Yb = cat_vol_morph(~Ycp | (Ycls{3}>128),'lc',1); - Ybd = cat_vbdist(single(~Yb),Yb,vx_vol); - Yvt = (Yg+abs(Ydiv))>0.4 & smooth3(single(Ycls{1})/255)<0.5 & Ybd>20 & ... + if debug==0, clear Yb; end + + + % Ysw .. save WM + % Ybv .. possible blood vessels + Ysw1 = cat_vol_morph(Ycls{2}>128 & (min(1,Ym)-Ydiv)<1.5,'lc',vxv*2) & (Ym - Ydiv)>5/6; + Ybv = ((min(1,Ym) - Ydiv + Yg)>2.0 | (Ycls{5}>16 & Ym<0.6 & Ycls{1}<192)) & ... + ~cat_vol_morph(Ysw1,'d',1) & Ym>0.2; + + + % Ycp = mask for CSF/BG distance initialization + Ycp1 = Ycls{3}>240 & Ydiv>0 & Yp0<1.1 & Ym<0.5 ; % typical CSF + Ycp2 = Ycls{5}>8 & Ycls{2}<32 & Ym<0.6 & Ydiv>0; % venes + Ycp3 = (Ym-Ydiv/4<0.4) & Ycls{3}>4 & Ycls{3}>16; % sulcal CSF + Ycp4 = (single(Ycls{6}) + single(Ycls{4}))>192; % non-CSF .. class 5 with error for negative t1 values + Ycp5 = ~cat_vol_morph(Ybb,'lc',5); % add background + Ycp6 = Ym<0.3; % but do not trust the brain mask! + Ycp = Ycp1 | Ycp2 | Ycp3 | Ycp4 | Ycp5 | Ycp6; % combine maps + Ycp(smooth3(Ycp)>0.4) = 1; % remove some meninges + if debug==0; clear Ycp1 Ycp2 Ycp3 Ycp4 Ycp5 Ycp6; end + + + % Ywd .. WM distance map + % Ysk .. divergence skeleton to improve CSF map + % deactivated due to problems in non-cortical structures + Ywd = cat_vbdist(single(Yp0>2.5),Yp0>0.5,vx_vol); % WM distance for skeleton + %Ysk = cat_vol_div(min(5,Ywd),2); % divergence skeleton + %Ysk = (Ym + min(0,Ysk))<0.2; % binary divergence skeleton + %Ycp = Ycp | Ysk; if debug==0, clear Ysk; end % combine skeleton with CSF map + + % Ycd .. CSF/BG distance map + Ycd = cat_vbdist(single(Ycp),~Ycp,vx_vol); % real CSF distance + Ycd((Ym-Ydiv<2/3 | Ydiv>0.1) & Ycls{3}>4 & Ycls{3}>1) = ... % correction for sulci + min(Ycd((Ym-Ydiv<2/3 | Ydiv>0.1) & Ycls{3}>4 & Ycls{3}>1),1.5); % ... maybe another distance estimation? + + + % we need to remove strong edge regions, because here is no GM layer between CSF and WM ??? + % Yb .. brain mask + % Ybd .. skull distance +%%% 20180801 - What is the difference of Yb2 to Ybb and the (previously replaced) Yb? + Yb2 = cat_vol_morph(~Ycp | (Ycls{3}>128),'lc',1); + Ybd = cat_vbdist(single(~Yb2),Yb2,vx_vol); + + + % Yvt .. Ventricle map WITHOUT atlas data as large deep CSF areas + % need the ventricle to avoid PVE GM between WM and CSF. + % There is an update for atlas data. + Yvt = (Yg + abs(Ydiv))>0.4 & smooth3(single(Ycls{1})/255)<0.5 & Ybd>20 & ... cat_vol_morph(Ycls{3}>8,'d',vxv) & cat_vol_morph(Ycls{2}>8,'d',vxv); Yvt = smooth3(Yvt)>0.7; Yvt = smooth3(Yvt)>0.2; + + %% final tissue maps: Ycm = CSF, Ygm = GM, Ywm = WM - Ysc = Ycp & Yb & Ycls{3}>192 & ~Ybv & Ym<0.45 & Yg<0.1; - Ycm = Ycp & Yb & Ycls{3}>192 & ~Ybv & (Yb | Ym>1/6) & Ym<0.45 & Yg<0.25 & Ym>0; % & Ydiv>-0.05; - if ~debug, clear Ycp; end - %Ycm = Ycm | (Yb & (Ym-max(0,Ydiv))<0.5); - Ywm = (Ysw | Ycls{2}>252 | ((Ycd-Ydiv)>2 & Ydiv<0 & Ym>0.9+LASstr*0.05 & Yb) | ... % save WM - ((Ycd-Ydiv.*Ycd)>4 & (Ydiv<-0.01) & Yb & Ym>0.5 & Ybd<20 & Ycd>2) ) & ... - ... ((Ycd-Ydiv*5)>3 & (Ydiv<-0.01 & (Yg + max(0,0.05-Ycd/100))<0.1) & Yb & Ym>0.4 & Ybd<20 & Ycd>2.5) ) & ... % further WM - ~Ybv & Yb & Ybd>1 & (Ycd>1.0 | (Yvt & Yp0>2.9)) & (Yg+Ydiv<(Ybd/50) | (Ydiv-Ym)<-1); % Ybd/800 + Ycd/50 - if ~debug, clear Ysw; end - Ygm = ~Yvt & Ybb & ~Ybv & ~Ywm & ~Ycm & Ycd>0.5 & (Ym-Ydiv-max(0,2-Ycd)/10)<0.9 & ... (Ym+Ydiv)>0.5 & ... ~Ysk & - (Ycls{1}>4 | (Ym>0.7 & Ycls{3}>64) | Ycd<(Ym+Ydiv)*3 ) & ... - smooth3(Yg>(Ybd/800) & Ycls{2}<240 )>0.6; % avoid GM next to hard boundies in the middle of the brain - if ~debug, clear Ybv Yvt; end - Ygx = Ybb & ~Ycm & ~Ywm & Ym>1/3 & Ym<2.8/3 & Yg<0.4 & (Ym-Ydiv)>1/3 & (Ym-Ydiv)<1; Ygx(smooth3(Ygx)<0.5) = 0; - Ygm = Ygm | Ygx; clear Ygx; - Ygm = Ygm | (Ym>1.5/3 & Ym<2.8/3 & ~Ywm & ~Ycm & Ybb); - Ygm(smooth3(Ygm)<0.25)=0; - %Ygw = Ygm & smooth3(Ywm)<0.1 & smooth3(Ycm)<0.4 & Ycd>0 & Ycd<2 & Ydiv<0.4 & Ydiv>-0.3 & Yg<0.1; %& (Ydiv>-0.4 | Ycd>1.5) - if ~debug, clear Ybv Ycp; end %Ycd + % ------------------------------------------------------------------- + + % CSF: + Ycm = Ycp & Yb2 & Ycls{3}>192 & ~Ybv & Ym<0.45 & Yg<0.25 & Ym>0 ; + %Ycm = Ycm | (Yb2 & (Ym - max(0,Ydiv))<0.5); % adding of divergence information was not save + if debug==0, clear Ycp; end - if debug>1 - try %#ok - [pth,nam] = spm_fileparts(res.image0(1).fname); tpmci=0; - tpmci=tpmci+1; tmpmat = fullfile(pth,reportfolder,sprintf('%s_%s%02d%s.mat',nam,'LAS',tpmci,'prepeaks')); - save(tmpmat); - end - end - %% + + %% WM: + % Different definitions of the possible WM (Ysw,Ysw2,Ysw3,Ysw4) + % and a variable (Ygw) to avoid non WM areas. + % Ysw1 .. defined above + % Ysw2 .. addition WM map that further include CSF distance (Ycd) and + % divergence information (Ydiv) and use a flexible intensity + % Ysw3 .. similar to Ysw2 with a skull-near and CSF distance criteria + % to reconstruct WM gyri + % Ysw4 .. was remove long ago + % Ygw .. map to of regions that we want to avoid in all possible WM + % definitions Ysw* + Ysw2 = Yb2 & (Ycd - Ydiv)>2 & Ydiv<0 & Ym>(0.9 + LASstr * 0.05); % general WM + Ysw3 = Yb2 & (Ycd - Ydiv .* Ycd)>4 & Ydiv<-0.01 & Ym>0.5 & Ybd<20 & Ycd>2; % addition WM gyri close to the cortex + %Ysw4 = ( (Ycd - Ydiv*5)>3 & Yb2 & Ym>0.4 & Ybd<20 & Ycd>2.5) & ... % further WM (rmoved long ago) + % (Ydiv<-0.01 & (Yg + max(0,0.05 - Ycd/100))<0.1 ); + + Ygw = ~Ybv & Yb2 & Ybd>1 & (Ycd>1.0 | (Yvt & Yp0>2.9)) & ... + (Yg + Ydiv<(Ybd/50) | (Ydiv - Ym)<-1); + + Ywm = ( Ycls{2}>252 | Ysw1 | Ysw2 | Ysw3 ) & Ygw; + if debug==0, clear Ysw Ysw2 Ysw3 Ygw; end + + %% GM: + % Different criteria (Ygm1,Ygm2,Ygm3) to define the possible GM area + % whereas Ygm4 describe general GM limitation. + % Ygm1 .. upper intensity limitation + % Ygm2 .. lower intensity limitation + % Ygm3 .. avoid GM next to hard boundaries in the middle of the brain + Ygm1 = ( Ym - Ydiv - max(0,2 - Ycd)/10 )<0.9; % upper limit + Ygm2 = Ycls{1}>4 | ( Ym>0.7 & Ycls{3}>64 ) | Ycd<(Ym + Ydiv)*3; % lower limit + Ygm3 = smooth3(Yg>(Ybd/800) & Ycls{2}<240 )>0.6; % avoid GM-WM PVE + Ygm = Ybb & ~Yvt & ~Ybv & ~Ywm & ~Ycm & Ycd>0.5 & Ygm1 & Ygm2 & Ygm3; + if debug==0, clear Ybv Yvt Ygm1 Ygm2 Ygm3; end + + % Ygm4 .. general GM limitations + Ygm4 = Ybb & ~Ycm & ~Ywm & Ym>1/3 & Ym<2.8/3 & ... + Yg<0.4 & (Ym - Ydiv)>1/3 & (Ym - Ydiv)<1; + Ygm4( smooth3(Ygm4)<0.5 ) = 0; % remove small dots + + % combine possible area with limitation map + Ygm = Ygm | Ygm4; + if debug==0, clear Ygm4; end + Ygm = Ygm | (Ym>1.5/3 & Ym<2.8/3 & ~Ywm & ~Ycm & Ybb); + Ygm(smooth3(Ygm)<0.25) = 0; + if debug==0, clear Ybv Ycp; end + + + + + + %% further maps that require the atlas if LABl1 - %% ------------------------------------------------------------------ % SPM GM segmentation can be affected by inhomogeneities and some GM - % is missclassified as CSF/GM (Ycls{5}). But for some regions we can - % trust these information more - % ------------------------------------------------------------------ - Ybd = cat_vbdist(single(~Yb),Yb,vx_vol); - Ycbp = cat_vbdist(single(NS(Yl1,LAB.CB)),Yb,vx_vol); % next to the cerebellum - Ycbn = cat_vbdist(single(~NS(Yl1,LAB.CB)),Yb,vx_vol); % not to deep in the cerebellum - Ylhp = cat_vbdist(single(mod(Yl1,2)==1 & Yb & Yl1>0),Yb,vx_vol); % GM next to the left hemisphere - Yrhp = cat_vbdist(single(mod(Yl1,2)==0 & Yb & Yl1>0),Yb,vx_vol); % GM next to the righ hemishpere - Ybv2 = Ycls{5}>2 & Ym<0.7 & Ym>0.3 & Yb & (... + % areas are miss-classified as CSF/GM (Ycls{5}) when the affine + % registration was not optimal. But for some regions we can trust + % these information. + + + + %% regions for cleanup of meninges and blood vessels + % Ycbp .. close to the cerebellum + % Ycbn .. not to deep in the cerebellum + % Ylhp / Yrhp .. GM next to the left/right hemisphere + % Ybv2 .. close to the skull and between the hemispheres or between cerebrum and cerebellum + % Ybvv .. sinus + Ycbp = cat_vbdist(single( NS(Yl1,LAB.CB)),Yb2,vx_vol); % close to cerebellum + Ycbn = cat_vbdist(single(~NS(Yl1,LAB.CB)),Yb2,vx_vol); % close to cerebellum surface + Ylhp = cat_vbdist(single(mod(Yl1,2)==1 & Yb2 & Yl1>0),Yb2,vx_vol); % GM next to the left hemisphere + Yrhp = cat_vbdist(single(mod(Yl1,2)==0 & Yb2 & Yl1>0),Yb2,vx_vol); % GM next to the right hemisphere + Ybv2 = Ycls{5}>2 & Ym<0.7 & Ym>0.3 & Yb2 & (... ((Ylhp+Ybd/2)0.5; - Ybvv = (Ym-max(0,6-abs(Ycbp-6))/50)<0.6 & Ym>0.4 & Yb & Ycbp<8 & Ycbp>1; - if ~debug, clear Ycbp; end + Ybvv = (Ym - max(0,6 - abs(Ycbp - 6))/50)<0.6 & Ym>0.4 & Yb2 & Ycbp<8 & Ycbp>1; % sinus + if debug==0, clear Ycbp; end + - %% subcortical map refinements + + %% refinements of subcortical regions + % Thalamus (TH): + % YTH .. enlarged atlas ROI + % Ytd .. CSF distance in the thalamus YTH + % Yxd .. distance to the basal-ganglia (BGL) in the thalamus YTH THth = 0.8 - LASstr*0.6; %0.5; % lower more thalamus - YTH = NS(Yl1,LAB.TH) | (cat_vol_morph(NS(Yl1,LAB.TH),'d',3) & Ym>0.5 & Ycls{1}>128); - Ytd = cat_vbdist(single(Ym<0.45),YTH | NS(Yl1,LAB.BG),vx_vol); Ytd(Ytd>2^16)=0; % CSF distance in the TH - Yxd = cat_vbdist(single(NS(Yl1,LAB.BG)),YTH,vx_vol); Yxd(Yxd>2^16)=0; % BG distance in the TH - %Yyd = cat_vbdist(single(NS(Yl1,LAB.TH)),NS(Yl1,LAB.BG),vx_vol); Yyd(Yyd>2^16)=0; % TH distance in the BG + YTH = NS(Yl1,LAB.TH) | (cat_vol_morph(NS(Yl1,LAB.TH),'d',3) & Ym>0.5 & Ycls{1}>128); + Ytd = cat_vbdist(single(Ym<0.45),YTH | NS(Yl1,LAB.BG),vx_vol); Ytd(Ytd>2^16)=0; % CSF distance in the TH + Yxd = cat_vbdist(single(NS(Yl1,LAB.BG)),YTH,vx_vol); Yxd(Yxd>2^16)=0; % BGL distance in the TH + if debug==0, clear YTH; end + + + % Yss .. (enlarged) subcortical structures + % Ynw .. no WM ?? Yss = NS(Yl1,LAB.BG) | NS(Yl1,LAB.TH); - Yss = Yss | (cat_vol_morph(Yss,'d',vxv*2) & Ym>2.25/3 & Ym<2.75/3 & Ydiv>-0.01); % add ihger tissue around mask - Yss = Yss | (cat_vol_morph(Yss,'d',vxv*3) & NS(Yl1,LAB.VT) & Yp0>1.5 & Yp0<2.3); % add lower tissue around mask - Yss = Yss & Yp0>1.5 & (Yp0<2.75 | (Ym<(2.5+LASstr*0.45)/3 & Ydiv>-0.05)); % by intensity - Yss = Yss | ((Yxd./max(eps,Ytd+Yxd))>THth/2 & (Yp0<2.75 | (Ym<(2.75+LASstr*0.20)/3 & Ydiv>-0.05))); % save TH by distances - for overcorrected images - Yss = cat_vol_morph(Yss,'o'); - Ynw = (Yxd./max(eps,Ytd+Yxd))>THth/2 | (NS(Yl1,LAB.BG) & Ydiv>-0.01); - if ~debug, clear Ytd Yxd ; end - % increase CSF roi + Yss = Yss | (cat_vol_morph(Yss,'d',vxv*2) & Ym>2.25/3 & Ym<2.75/3 & Ydiv>-0.01); % add higher tissue around mask + Yss = Yss | (cat_vol_morph(Yss,'d',vxv*3) & NS(Yl1,LAB.VT) & Yp0>1.5 & Yp0<2.3); % add lower tissue around mask + Yss = Yss & Yp0>1.5 & (Yp0<2.75 | (Ym<(2.5 + LASstr * 0.45)/3 & Ydiv>-0.05)); % limit the enlarged map by intensity + Yss = Yss | ((Yxd./max(eps,Ytd + Yxd))>THth/2 & ... % save TH by distances - for overcorrected images + (Yp0<2.75 | (Ym<(2.75 + LASstr * 0.20)/3 & Ydiv>-0.05))); + Yss = cat_vol_morph(Yss,'o'); + + Ynw = (Yxd./max(eps,Ytd + Yxd))>THth/2 | (NS(Yl1,LAB.BG) & Ydiv>-0.01); + if debug==0, clear Ytd Yxd ; end + + + % increase CSF ROI + % Yvt .. improve ventricle ROI with atlas information +%%%%% 20180801 - Why not use the ventricle maps from the partitioning? + % - Because the partitioning is after LAS + % - Why not use a fast partitioning before LAS? + % Why not use the code from the partitioning? + % Is a good ventricle segmentation here not important? Yvt = cat_vol_morph( (NS(Yl1,LAB.VT) | cat_vol_morph(Ycm,'o',3) ) ... - & Ycm & ~NS(Yl1,LAB.BG) & ~NS(Yl1,LAB.TH) & Ybd>30,'d',vxv*3) & ~Yss; % ventricle roi to avoid PVE GM between WM and CSF - Ycx = (NS(Yl1,LAB.CB) & ((Ym-Ydiv)<0.55 | Ycls{3}>128)) | (((Ym-Ydiv)<0.45 & Ycls{3}>8)| Ycls{3}>240); - % in the crebellum tissue can be differentated by div etc. - Ycwm = NS(Yl1,LAB.CB) & (Ym-Ydiv*4)>5/6 & Ycd>3 & Yg>0.05; + & Ycm & ~NS(Yl1,LAB.BG) & ~NS(Yl1,LAB.TH) & Ybd>30,'d',vxv*3) & ~Yss; + + + + %% meningeal corrections + % Ycx .. cerebellar and other CSF + % Ycwm .. cerebellar WM + % Yccm .. cerebellar CSF + % Ybwm .. brain WM + % Ybcm .. brain CSF + Ycx = ( NS(Yl1,LAB.CB) & ( (Ym - Ydiv)<0.55 | Ycls{3}>128 ) ) | ... + ( ( (Ym - Ydiv)<0.45 & Ycls{3}>8 ) | Ycls{3}>240 ); + % in the cerebellum tissue can be differentiated by div etc. + Ycwm = NS(Yl1,LAB.CB) & (Ym - Ydiv*4)>5/6 & Ycd>3 & Yg>0.05; Yccm = NS(Yl1,LAB.CB) & Ydiv>0.02 & Ym<1/2 & Yg>0.05; - Ybwm = (Ym-Ydiv*4)>0.9 & Ycd>3 & Yg>0.05; %Ydiv<-0.04 & Ym>0.75 & Ycd>3; + Ybwm = (Ym - Ydiv*4)>0.9 & Ycd>3 & Yg>0.05; Ybcm = Ydiv>0.04 & Ym<0.55 & Yg>0.05; + + + %% correction 1 of tissue maps - Ywmtpm = ( single( Ywtpm )/255 .*Ym.*(1-Yg-Ydiv).*cat_vol_morph(NS(Yl1,1).*Ybd/5,'e',1))>0.6; % no WM hyperintensities in GM! - if ~debug, clear Ywtpm; end - % + % Ywmtpm .. refined WM template map (original Ywtpm) + Ywmtpm = ( single( Ywtpm )/255 .* Ym .* (1 - Yg - Ydiv) .* ... + cat_vol_morph( NS(Yl1,1) .* Ybd/5 ,'e',1) )>0.6; % no WM hyperintensities in GM! + if debug==0, clear Ywtpm; end + + + % refine the GM map Ygm = Ygm | (Yss & ~Yvt & ~Ycx & ~Ybv2 & ~Ycwm & ~(Yccm | Ybcm)); Ygm = Ygm & ~Ywmtpm & ~Ybvv; % no WMH area - Ywm = (Ywm & ~Yss & ~Ybv2 & ~Ynw) | Ycwm | Ybwm; %& ~NS(Yl1,LAB.BG) - Ywmtpm(smooth3(Ywmtpm & Ym<11/12)<0.5)=0; + + + % refine the WM map + Ywm = Ywm & ~Yss & ~Ybv2 & ~Ynw; + Ywm = Ywm | Ycwm | Ybwm; + Ywmtpm(smooth3(Ywmtpm & Ym<11/12)<0.5) = 0; % 20180801 - Why here? Ywm = Ywm & ~Ywmtpm & ~Ybvv & ~Yss; % no WM area - Ycm = Ycm | ( (Ycx | Yccm | Ybcm) & Yg<0.2 & Ym>0 & Ydiv>-0.05 & Ym<0.3 & Yb ) | Ybvv; - if ~debug, clear Ycwm Yccm Ycd Ybvv; end - % mapping of the brainstem to the WM (well there were some small GM + + + Ycm = Ycm | ( (Ycx | Yccm | Ybcm) & Yg<0.2 & Ym>0 & Ydiv>-0.05 & Ym<0.3 & Yb2 ) | Ybvv; + if debug==0, clear Ycwm Yccm Ycd Ybvv Ycx; end + + + % Mapping of the brain stem to the WM (well there were some small GM % structures, but the should not effect the local segmentation to much. - Ybs = cat_vol_morph(NS(Yl1,LAB.BS) & Ym<1.2 & Ym>0.9 & Yp0>2.5,'c',2*vxv) & Ym<1.2 & Ym>0.9 & Yp0>1.5; + % Ybs .. brain stem mask + Ybs = Ym<1.2 & Ym>0.9 & Yp0>1.5 & ... + cat_vol_morph(NS(Yl1,LAB.BS) & Ym<1.2 & Ym>0.9 & Yp0>2.5,'c',2*vxv); Ygm = (Ygm & ~Ybs & ~Ybv2 & ~Ywm) | Yss; Ywm = Ywm | (Ybs & Ym<1.1 & Ym>0.9 & Yp0>1.5) ; - if ~debug, clear Ycx; end end - % back to original resolution for full bias field estimation + + + + + + + %% back to original resolution for full bias field estimation + % ------------------------------------------------------------------- Ycls = Yclso; clear Yclso; Ycm = cat_vol_resize(Ycm , 'dereduceBrain' , BB); Ygm = cat_vol_resize(Ygm , 'dereduceBrain' , BB); Ywm = cat_vol_resize(Ywm , 'dereduceBrain' , BB); Yvt = cat_vol_resize(Yvt , 'dereduceBrain' , BB); - Yb = cat_vol_resize(Yb , 'dereduceBrain' , BB); + Yb2 = cat_vol_resize(Yb2 , 'dereduceBrain' , BB); Yss = cat_vol_resize(Yss , 'dereduceBrain' , BB); Ybb = cat_vol_resize(Ybb , 'dereduceBrain' , BB); - Ysc = cat_vol_resize(Ysc , 'dereduceBrain' , BB); Ybs = cat_vol_resize(Ybs , 'dereduceBrain' , BB); Ybv2 = cat_vol_resize(Ybv2 , 'dereduceBrain' , BB); @@ -325,256 +509,394 @@ Yg = cat_vol_resize(Yg , 'dereduceBrain' , BB); Ydiv = cat_vol_resize(Ydiv , 'dereduceBrain' , BB); Ywd = cat_vol_resize(Ywd , 'dereduceBrain' , BB); + - if debug>1 - try %#ok % windows requires this... i don't know why ... maybe the file size - tpmci=tpmci+1; tmpmat = fullfile(pth,reportfolder,sprintf('%s_%s%02d%s.mat',nam,'LAS',tpmci,'prepeaks')); - save(tmpmat); - end - end - % correction for negative values - [Ysrcx,thx] = cat_stat_histth(Ysrc,99); clear Ysrcx; - srcmin = thx(1); %min(Ysrc(:)); - Ysrc = Ysrc - srcmin; T3th = T3th - srcmin; Tthc = Tth; Tthc.T3th = Tth.T3th - srcmin; + [Ysrcx,thx] = cat_stat_histth(Ysrc,99); clear Ysrcx; %#ok + srcmin = thx(1); + Ysrc = Ysrc - srcmin; + T3th = T3th - srcmin; + Tthc = Tth; Tthc.T3th = Tth.T3th - srcmin; if exist('Ysrco2','var'),Ysrco2 = Ysrco2 - srcmin; end + + + + + + %% --------------------------------------------------------------------- % Now, we can estimate the local peaks % --------------------------------------------------------------------- - % Estimation of the local WM threshold with "corrected" GM voxels to - % avoid overfitting (see BWP cerebellum). - % CSF is problematic in high contrast or skull-stripped image should - % not be used here, or in GM peak estimation - mres = 1.1; - stime = cat_io_cmd(' Estimate local tissue thresholds (WM)','g5','',verb,stime); dispc=dispc+1; +% Estimation of the local WM intensity with help of intensity corrected +% GM, CSF and head tissues to avoid overfitting (eg. BWP cerebellum). +% CSF is problematic in high contrast or skull-stripped images and +% should not be used here or in the GM peak estimation. +% --------------------------------------------------------------------- + mres = 1.1; + stime = cat_io_cmd(' Estimate local tissue thresholds (WM)','g5','',verb,stime); Ysrcm = cat_vol_median3(Ysrc.*Ywm,Ywm,Ywm); rf = [10^5 10^4]; - T3th3 = max(1,min(10^6,rf(2) / max(eps,round(T3th(3)*rf(1))/rf(1)))); + T3th3 = max(1,min(10^6,rf(2) / (round(T3th(3)*rf(1))/rf(1)))); Ysrcm = round(Ysrcm*T3th3)/T3th3; - Ygw2 = Ycls{1}>128 & Ym>2/3-0.04 & Ym<2/3+0.04 & Ygm .*Ydiv>0.01; - Ygw2 = Ygw2 | (Ycls{1}>128 & Yg<0.05 & abs(Ydiv)<0.05 & ~Ywm & Ym<3/4); % large stable GM areas - like the BWP cerebellum - Ygw3 = Ycls{3}>128 & Yg<0.05 & ~Ywm & ~Ygm & Ywd<3; - Ygw3(smooth3(Ygw3)<0.5)=0; - [Yi,resT2] = cat_vol_resize(Ysrcm,'reduceV',vx_vol,mres,32,'max'); % maximum reduction for the WM - % - if cat_stat_nanmean(Ym(Ygw3))>0.1, % not in images with to low CSF intensity (error in skull-stripped) - Ygi = cat_vol_resize(Ysrc.*Ygw2*T3th(3)/mean(Ysrc(Ygw2(:))) + Ysrc.*Ygw3*T3th(3)/mean(Ysrc(Ygw3(:))),'reduceV',vx_vol,mres,32,'meanm'); % mean for other tissues - else - Ygi = cat_vol_resize(Ysrc.*Ygw2*T3th(3)/mean(Ysrc(Ygw2(:))),'reduceV',vx_vol,mres,32,'meanm'); % mean for other tissues + + + % Ygwg .. large homogen GM regions (e.g. subcortical structures or in the BWP cerebellum) + % Ygwc .. large homogen CSF regions + Ygwg = Ycls{1}>128 & Ym>(2/3 - 0.04) & Ym<(2/3 + 0.04) & Ygm .*Ydiv>0.01; + Ygwg = Ygwg | (Ycls{1}>128 & Yg<0.05 & abs(Ydiv)<0.05 & ~Ywm & Ym<3/4); + Ygwc = Ycls{3}>128 & Yg<0.05 & ~Ywm & ~Ygm & Ywd<3; + Ygwc(smooth3(Ygwc)<0.5) = 0; + if debug==0, clear Ywd; end + + + % also use normalized GM tissues + % Ygi .. tissue-corrected intensity map (normalized for WM) + Ygi = Ysrc .* Ygwg * T3th(3)/mean(Ysrc(Ygwg(:))); + if cat_stat_nanmean(Ym(Ygwc))>0.1 + % use normalized CSF tissue only in images with clear CSF intensity + % (to avoid errors in skull-stripped data) + Ygi = Ygi + Ysrc .* Ygwc * T3th(3)/mean(Ysrc(Ygwc(:))); end - for xi=1:2*LASi, Ygi = cat_vol_localstat(Ygi,Ygi>0,2,1); end; Ygi(smooth3(Ygi>0)<0.3)=0; - Yi = cat_vol_localstat(Yi,Yi>0,1,3); % one maximum for stabilization of small WM structures - Yi(Yi==0 & Ygi>0)=Ygi(Yi==0 & Ygi>0); - for xi=1:2*LASi, Yi = cat_vol_localstat(Yi,Yi>0,2,1); end % no maximum here! - % - if cat_stat_nanmean(Ym(Ygw3))>0.1 && cat_stat_nanmean(Ysrc(Ycls{6}(:)>128))>T3th(1) - %% - Ygw4 = Ycls{6}>240 & ~Ygw3 & ~Ygw2 & Yg<0.5 & abs(Ydiv)<0.1 & ... - Ysrc>min(res.mn(res.lkp==6))*0.5 & Ysrc0,2,1); end; % intensity smoothing in the GM area + Ygi(smooth3(Ygi>0)<0.3) = 0; % remove small dots (meninges) + Yi = cat_vol_localstat(Yi,Yi>0,1,3); % maximum filtering to stabilize small WM structures + Yi(Yi==0 & Ygi>0) = Ygi(Yi==0 & Ygi>0); % combine the WM and GM map + for xi = 1:2*LASi, Yi = cat_vol_localstat(Yi,Yi>0,2,1); end % intensity smoothing in both regions no maximum here! + + + % Add head tissue if it is available. + if cat_stat_nanmean(Ym(Ygwc))>0.1 && cat_stat_nanmean(Ysrc(Ycls{6}(:)>128))>T3th(1) + try % no skull-stripped + % definion of head tissue + Ygwh = Ycls{6}>128 & ~Ygwc & ~Ygwg & Yg<0.5 & abs(Ydiv)<0.1 & ... + Ysrc>min( res.mn(res.lkp==6) * 0.5 ) & Ysrc ( cat_stat_nanmedian( Ysrc(Ygwh(:)) ) + 2 * cat_stat_nanstd( Ysrc(Ygwh(:)) ) ); + Ygwh( Ygwhn ) = false; clear Ygwhn + %% go to low resolution + [Ygi,resTB] = cat_vol_resize(Ysrc .* Ygwh * T3th(3)/mean(Ysrc(Ygwh(:))),'reduceV',vx_vol,mres,32,'meanm'); + % additional approximation + Ygi = cat_vol_approx(Ygi,'nh',resTB.vx_volr,2); Ygi = cat_vol_smooth3X(Ygi,2 * LASfs); + Ygi = Ygwh .* max(eps,cat_vol_resize(Ygi,'dereduceV',resTB)); + Yi(Yi==0) = Ygi(Yi==0); + if debug==0; clear Ygwh; end + end end - if debug==0; clear Ygw2 Ygw3 Ygw4; end - Yi(Yi==0) = Ygi(Yi==0); - Yi = cat_vol_approx(Yi,'nh',resT2.vx_volr,2); Yi = cat_vol_smooth3X(Yi,2*LASfs); + if debug==0, clear Ygwg Ygwc; end + + + % final approximation of the WM inhomogeneity + Yi = cat_vol_approx(Yi,'nh',resT2.vx_volr,2); + Yi = cat_vol_smooth3X(Yi,2 * LASfs); Ylab{2} = max(eps,cat_vol_resize(Yi,'dereduceV',resT2)); - % Ylab{2} = Ylab{2} .* mean( [median(Ysrc(Ysw(:))./Ylab{2}(Ysw(:))),1] ); +% Ylab{2} = Ylab{2} .* mean( [median(Ysrc(Ysw(:))./Ylab{2}(Ysw(:))),1] ); if debug==0; clear Ysw Yi; end - %% update GM tissue map - %stime = cat_io_cmd(' Update tissues maps','g5','',verb,stime); dispc=dispc+1; - %Ybb = cat_vol_morph((Yb & Ysrc./Ylab{2}<(T3th(1) + 0.25*diff(T3th(1:2))) & Ydiv<0.05) | Yp0>1.5,'lo',1); - Ygm(Ysrc./Ylab{2}>(T3th(2) + 0.90*diff(T3th(2:3)))/T3th(3))=0; % correct GM mean(T3th([2:3,3]))/T3th(3) - Ygm(Ysrc./Ylab{2}<(T3th(2) + 0.75*diff(T3th(2:3)))/T3th(3) & ... - Ysrc./Ylab{2}<(T3th(2) - 0.75*diff(T3th(2:3)))/T3th(3) & ... - Ydiv<0.3 & Ydiv>-0.3 & Ybb & ~Ywm & ~Yvt & ~Ybv2 & Ycls{1}>48)=1; - Ywmd2 = cat_vbdist(single(Ywm),Yb); - Ygx = Ywmd2-Ym+Ydiv>0.5 & Ym+0.5-Ydiv-Yg-Ywmd2/10>1/3 & ~Ybv2 & ... low intensity tissue - ~(Ym-min(0.2,Yg+Ywmd2/10-Ydiv)<1/4) & Yg0.5 & ~Ywm; - Ygx(smooth3(Ygx)<0.5)=0; - Ygm = Ygm | Ygx; % correct gm (spm based) - if debug==0; clear Ygx; end - - %% - Ycm = ~Ygm & ~Ywm & ~Ybv2 & Yg<0.6 & (Ycm | (Yb & (Ysrc./Ylab{2})<((T3th(1)*0.5 + 0.5*T3th(2))/T3th(3)))); - % - Ycp = (Ycls{2}<128 & Ydiv>0 & Yp0<2.1 & Ysrc./Ylab{2}32 & Ycls{2}<32 & Ysrc./Ylab{2}0) | ... % venes - ((Ym-Ydiv<0.4) & Ycls{3}>4 & Ycls{3}>16 & Ysrc./Ylab{2}192 | ... % save non-csf - Ysrc./Ylab{2}0.4)=1; % remove some meninges - Ycd = cat_vbdist(single(Ycp),~Ycp,vx_vol); - if debug==0; clear Ycp; end - %% - Ygm = Ygm & ~Ycm & ~Ywm & Ywd<5; % & ~Ybvv & ~Ysk - Ygm = Ygm | (NS(Yl1,1) & Ybd<20 & (Ycd-Ydiv)<2 & Ycls{1}>0 & ~Ycm & Ybb & Ym>0.6 & Yg(T3th(2) + 0.90 * diff(T3th(2:3)))/T3th(3)) = 0; + Ygm(Ysrc./Ylab{2}<(T3th(2) + 0.75 * diff(T3th(2:3)))/T3th(3) & ... + Ysrc./Ylab{2}<(T3th(2) - 0.75 * diff(T3th(2:3)))/T3th(3) & ... + Ydiv<0.3 & Ydiv>-0.3 & Ybb & ~Ywm & ~Yvt & ~Ybv2 & Ycls{1}>48) = 1; + Ywd2 = cat_vbdist(single(Ywm),Yb2); % update WM distance + Ygmn = (Ywd2 - Ym + Ydiv)>0.5 & ~Ybv2 & ... + (Ym + 0.5 - Ydiv - Yg - Ywd2/10)>1/3 & ... low intensity tissue + ~(Ym - min(0.2,Yg + Ywd2/10 - Ydiv)<1/4) & ... + Yg0.5 & ~Ywm; + Ygmn(smooth3(Ygmn)<0.5) = 0; + Ygm = Ygm | Ygmn; % correct gm (SPM based) + if debug==0; clear Ygm4; end + + + + + + + %% update maps +%%% ... more details + + % Update CSF maps. + % Ycm .. updated CSF map + % Ycp1 .. typical CSF + % Ycp2 .. venes + % Ycp3 .. sulcal CSF + % Ycp4 .. save non-CSF + % Ycp5 .. but do not trust the brain mask! + % Ycpn .. updated CSF 2 map + % Ycd2 .. updated CSF distance + Ycm = Ycm | (Yb2 & (Ysrc./Ylab{2})<((T3th(1)*0.5 + 0.5*T3th(2))/T3th(3))); + Ycm = ~Ygm & ~Ywm & ~Ybv2 & Yg<0.6 & Ycm; + Ycp1 = Ycls{2}<128 & Ydiv>0 & Yp0<2.1 & Ysrc./Ylab{2}32 & Ycls{2}<32 & Ysrc./Ylab{2}0; % venes + Ycp3 = (Ym - Ydiv<0.4) & Ycls{3}>4 & Ycls{3}>16 & Ysrc./Ylab{2}192; % save non-CSF + Ycp5 = Ysrc./Ylab{2}0.4 ) = 1; % remove some meninges + Ycd2 = cat_vbdist(single(Ycpn),~Ycpn,vx_vol); + + + %% update GM map + Ygm = Ygm & ~Ycm & ~Ywm & Ywd2<5; + if debug==0; clear Ywd2; end + Ygm = Ygm | (NS(Yl1,1) & Ybd<20 & (Ycd2 - Ydiv)<2 & ... + Ycls{1}>0 & ~Ycm & Ybb & Ym>0.6 & Yg1.5 | (Ym>1.2/3 & Ym<3.1/3 & Yb))>0.6,'lo',min(1,vxv)); % clear Yp0 Yvt - Ygm(~Ybb)=0; Ygm(smooth3(Ygm)<0.3)=0; - Ygm(smooth3(Ygm)>0.4 & Ysrc./Ylab{2}>mean(T3th(1)/T3th(2)) & Ysrc./Ylab{2}<(T3th(2)*0.2+0.8*T3th(3)))=1; - if debug==0; clear Ybb Yl1; end %Ydiv Yg + Ybb = cat_vol_morph(smooth3(Ygm | Ywm | Yp0>1.5 | ... + (Ym>1.2/3 & Ym<3.1/3 & Yb2))>0.6,'lo',min(1,vxv)); + if debug==0, clear Yp0 Yvt; end + Ygm(~Ybb) = 0; + Ygm(smooth3(Ygm)<0.3) = 0; + Ygm(smooth3(Ygm)>0.4 & Ysrc./Ylab{2}>mean(T3th(1)/T3th(2)) & ... + Ysrc./Ylab{2}<(T3th(2)*0.2+0.8*T3th(3))) = 1; + if debug==0; clear Ybb Yl1; end - stime = cat_io_cmd(' Estimate local tissue thresholds (GM)','g5','',verb,stime); dispc=dispc+1; + + + %% GM -% Yi = (Ysrc./Ylab{2} .* Ygm , Ysrc./Ylab{2}.*Ywm.*Ybs.*(T3th(2)+0.8*diff(T3th(2:3)))/T3th(3)); -% Yi(Ybv2) = Ysrc(Ybv2)./Ylab{2}(Ybv2) .* T3th(2)/mean(T3th(1:2)); + % .. more details + + % Yi .. main GM intensity map + stime = cat_io_cmd(' Estimate local tissue thresholds (GM)','g5','',verb,stime); Yi = Ysrc./Ylab{2} .* Ygm; - Yi = round(Yi*rf(2))/rf(2); - % Yi(Ybv2) = Ysrc(Ybv2)./Ylab{2}(Ybv2) .* T3th(2)/mean(T3th(1:2)); % ???? - Yi(Ybs) = Ysrc(Ybs)./Ylab{2}(Ybs) .* T3th(2)/T3th(3); + Yi = round(Yi * rf(2))/rf(2); + Yi(Ybs) = Ysrc(Ybs)./Ylab{2}(Ybs) .* T3th(2)/T3th(3); if debug==0; clear Ybs; end - Yi = cat_vol_median3(Yi,Yi>0.5,Yi>0.5); - Ycmx = smooth3(Ycm & Ysrc<(T3th(1)*0.8+T3th(2)*0.2))>0.9; Tcmx = mean(Ysrc(Ycmx(:))./Ylab{2}(Ycmx(:)))*T3th(3); + Yi = cat_vol_median3(Yi,Yi>0.5,Yi>0.5); % reduce artifacts + + % add CSF and CSF/GM areas (venes) to avoid overfitting + % Ycmx .. venes / sinus + % Tcmx .. intensity correction + Ycmx = smooth3( Ycm & Ysrc<(T3th(1)*0.8 + T3th(2)*0.2) )>0.9; + Tcmx = mean( Ysrc(Ycmx(:)) ./ Ylab{2}(Ycmx(:)) ) * T3th(3); Yi(Ycmx) = Ysrc(Ycmx)./Ylab{2}(Ycmx) .* T3th(2)/Tcmx; if debug==0; clear Ycm Ycmx; end - %Yii = Ysrc./Ylab{2} .* Ycm * T3th(2) / cat_stat_nanmedian(Ysrc(Ycm(:))); - [Yi,Yii,Ybgx,resT2] = cat_vol_resize({Yi,Ylab{2}/T3th(3),Ycls{6}>240},'reduceV',vx_vol,1,32,'meanm'); - for xi=1:2*LASi, Yi = cat_vol_localstat(Yi,Yi>0,3,1); end - Yi = cat_vol_approx(Yi,'nh',resT2.vx_volr,2); - Yi = min(Yi,Yii*(T3th(2) + 0.90*diff(T3th(2:3)))/T3th(3)); - Yi(Ybgx) = Yii(Ybgx)*cat_stat_nanmean(Yi(~Ybgx(:))); + + %% reduce data to lower resolution + Ybgx = Ycls{6}>128; + Ybgxn = Ysrc < ( cat_stat_nanmedian( Ysrc(Ybgx(:)) ) - 2 * cat_stat_nanstd( Ysrc(Ybgx(:)) ) ) | ... + Ysrc > ( cat_stat_nanmedian( Ysrc(Ybgx(:)) ) + 2 * cat_stat_nanstd( Ysrc(Ybgx(:)) ) ); + Ybgx = Ybgx & ~Ybgxn; clear Ybgxn; + [Yi,resT2] = cat_vol_resize( Yi ,'reduceV',vx_vol,1,32,'meanm'); + Yii = cat_vol_resize( Ylab{2}/T3th(3) ,'reduceV',vx_vol,1,32,'meanm'); + Ybgx = cat_vol_resize( Ybgx ,'reduceV',vx_vol,1,32,'meanm'); + for xi = 1:2*LASi, Yi = cat_vol_localstat(Yi,Yi>0,3,1); end % local smoothing + Yi = cat_vol_approx(Yi,'nh',resT2.vx_volr,2); % approximation of the GM bias + Yi = min(Yi,Yii*(T3th(2) + 0.90*diff(T3th(2:3)))/T3th(3)); % limit upper values to avoid wholes in the WM + Yi(Ybgx) = Yii(Ybgx) * cat_stat_nanmean(Yi(~Ybgx(:))); % use WM bias field in the background if debug==0; clear Yii Ybgx; end - %% - Yi = cat_vol_smooth3X(Yi,LASfs); - Ylab{1} = cat_vol_resize(Yi,'dereduceV',resT2).*Ylab{2}; + Yi = cat_vol_smooth3X(Yi,LASfs); % final low res smoothing + Ylab{1} = cat_vol_resize(Yi,'dereduceV',resT2).*Ylab{2}; % back to original resolution if debug==0; clear Yi; end - %Ylab{1}(Ygm) = Ysrc(Ygm); Ylab{1} = cat_vol_smooth3X(Ylab{1},LASfs); % can lead to overfitting - Ycm = (single(Ycls{3})/255 - Yg*4 + abs(Ydiv)*2)>0.5 & Ysrc<(Ylab{1}*mean(T3th([1,1:2]))/T3th(2)); %Ycm & Ysrc<(Ylab{1}*mean(T3th(1:2))/T3th(2)) & Yg<0.1; + + + % update CSF map + Ycm = (single(Ycls{3})/255 - Yg*4 + abs(Ydiv)*2)>0.5 & ... + Ysrc<(Ylab{1}*mean(T3th([1,1:2]))/T3th(2)); if debug==0; clear Ydiv; end Ycm(smooth3(Ycm)<0.5)=0; - Ycm(Yb & cat_vol_morph(Ysrc128 | (~Yb & Yg<=0.001),'e',4*vxv); - else - Ynb = cat_vol_morph(smooth3(Ycls{6})>128 | (~Yb & Yg<=0.001),'e',4*vxv); + % The differentiation of CSF and Background is not allway good. In some + % images with low intensity (especially for defacing/skull-stripping) + % the CSF can be similar to the background. + stime = cat_io_cmd(' Estimate local tissue thresholds (CSF/BG)','g5','',verb,stime); + Ynb = smooth3( Ycls{6})>128 & Yg<0.15; % & Ysrc0.5) * rf(2))/rf(2),'reduceV',vx_vol,8,16,'min');% only pure CSF !!! - [Yx,Yxm] = cat_vol_resize({round(Ysrc./Ylab{2} .* Ynb * rf(2))/rf(2),single(Ynb)},'reduceV',vx_vol,8,16,'mean');% only pure CSF !!! - Yx(Yxm<1) = 0; - %% + if debug==0, clear Ycpn; end + Yc = round(Ysrc./Ylab{2} .* (smooth3(Ycm)>0.5) * rf(2))/rf(2); + Yx = round(Ysrc./Ylab{2} .* Ynb * rf(2))/rf(2); + + % The correction should be very smooth (and fast) for CSF and we can use + % much lower resolution. + [Yc,resT2] = cat_vol_resize(Yc,'reduceV',vx_vol,8,16,'min'); % only pure CSF !!! + Yx = cat_vol_resize(Yx,'reduceV',vx_vol,8,16,'meanm'); + + % correct Nans and remove outlier + Yc(isnan(Yc)) = 0; Yc = cat_vol_median3(Yc,Yc~=0,true(size(Yc)),0.1); + Yx(isnan(Yx)) = 0; Yx = cat_vol_median3(Yx,Yx~=0,true(size(Yc)),0.1); + + % CSF intensity show be higher than background intensity Yx(Yc>0)=0; Yc(Yx>0)=0; - if cat_stat_nanmean(Yx(Yx~=0)) < cat_stat_nanmean(Yc(Yc~=0)) - meanYx = min(median(Yc(Yc(:)>0)),median(Yx(Yx(:)>0))); - meanYc = max(median(Yc(Yc(:)>0)),median(Yx(Yx(:)>0))); + if 0 % no so important ... keep it simple? + if cat_stat_nanmean(Yx(Yx~=0)) < cat_stat_nanmean(Yc(Yc~=0)) + meanYx = min(median(Yc(Yc(:)>0)),median(Yx(Yx(:)>0))); + meanYc = max(median(Yc(Yc(:)>0)),median(Yx(Yx(:)>0))); + else + meanYc = min(median(Yc(Yc(:)>0)),median(Yx(Yx(:)>0))); + meanYx = max(median(Yc(Yc(:)>0)),median(Yx(Yx(:)>0))); + end else - meanYc = min(median(Yc(Yc(:)>0)),median(Yx(Yx(:)>0))); - meanYx = max(median(Yc(Yc(:)>0)),median(Yx(Yx(:)>0))); - end - stdYbc = mean([std(Yc(Yc(:)>0)),std(Yx(Yx(:)>0))]); - %Yx = min(max(meanYx/2,Yx),min(meanYx*4,meanYc/2)); - %Yc = min(max(meanYc/2,Yx),meanYc/2); - Yxa = cat_vol_approx(Yx ,'nh',resT2.vx_volr,16); %+(Yb>0).*stdYbc + Yc.*meanYb/max(eps,meanYc) - Yca = cat_vol_approx(Yc + min(max( meanYx + stdYbc , meanYc - stdYbc ),... - Yx.*meanYc/max(eps,meanYx)),'nh',resT2.vx_volr,16); % + Yb.*meanYc/max(eps,meanYb) - if ~debug, clear Yc Yx; end - Yca = Yca*0.7 + 0.3*max(mean(Yca(:)),T3th(1)/T3th(3)); - %% - Yxa = cat_vol_smooth3X(Yxa,LASfs*2); - Yca = cat_vol_smooth3X(Yca,LASfs*2); - Ylab{3} = cat_vol_smooth3X(cat_vol_resize(Yca,'dereduceV',resT2).*Ylab{2},LASfs*2); - Ylab{6} = cat_vol_smooth3X(cat_vol_resize(Yxa,'dereduceV',resT2).*Ylab{2},LASfs*2); - if ~debug, clear Yxa Yca; end + meanYx = cat_stat_nanmean( [ 0 ; Yx(Yx(:)>0) ]); + meanYc = cat_stat_nanmedian( [ T3th(1) ; Yc(Yc(:)>0) ]); + end + % correct for high intensity backgrounds (e.g. Gabi's R1) + %Yx = Yx .* (Yx > (meanYc + max(0.1,std(Yc(Yc(:)>0))))); + + %% avoid CSF in background and background in CSF + if ~res.ppe.affreg.highBG + bcd = max( abs( [ diff( [meanYx meanYc] ) min( diff( T3th/max(T3th(:)) ) ) ] )); + Yxc = Yx + (Yx==0 & Yc~=0) .* (Yc * meanYx/meanYc - 0.5 * bcd) + ... % avoid BG in the CSF (white dots in the blue CSF) + ( (Yx~=0) * 0.2 * bcd ); % this term is the lower HG intensity (some nice noise in the BG) + if 0% (meanYx/meanYc) > 0.8 % this looks nice but its no good and change thickness + Ycc = Yc + (Yc==0 & Yx~=0) .* (Yx * meanYc/meanYx + 0.5 * bcd); % here we avoid CSF in the BG (higher boudnary) ... blue CSF dots in the BG + else + Ycc = Yc; + end + %Ycc = Yc + min(max( meanYx + stdYbc , meanYc - stdYbc ), Yx .* meanYc/max(eps,meanYx)); + else + Yxc = Yx; + Ycc = Yc; + end + %% guaranty small difference between backgound and CSF intensity + Yxa = cat_vol_approx(Yxc ,'nh',resT2.vx_volr,16); + Yca = cat_vol_approx(Ycc ,'nh',resT2.vx_volr,16); + if debug==0, clear Yc Yx; end + Yca = Yca * 0.7 + 0.3 * max( mean(Yca(:)) , T3th(1)/T3th(3) ); + % smoothing + Yxa = cat_vol_smooth3X( Yxa , LASfs * 2 ); + Yca = cat_vol_smooth3X( Yca , LASfs * 2 ); + %% back to original resolution, combination with WM bias, and final smoothing + Yca = cat_vol_resize(Yca,'dereduceV',resT2) .* Ylab{2}; + Yxa = cat_vol_resize(Yxa,'dereduceV',resT2) .* Ylab{2}; + Ylab{3} = cat_vol_smooth3X( Yca , LASfs * 2 ); + Ylab{6} = cat_vol_smooth3X( Yxa , LASfs * 2 ); + if debug==0, clear Yxa Yca; end + + if res.ppe.affreg.highBG + Ylab{6} = min(Ysrc(:)); + end + else + % simple + Ynb = smooth3( Ycls{6})>128 & Yg<0.3; + %Ynb = cat_vol_morph(Ynb,'e',4*vxv); + + Ylab{3} = ones(size(Ygm),'single') .* min( [ T3th(1) , cat_stat_nanmean(Ysrc(Ycm(:))) ] ); + Ylab{6} = ones(size(Ynb),'single') .* min( [ T3th(1) - min(diff(T3th)) , cat_stat_nanmean(Ysrc(Ycm(:))) ] ); + end + + % final scaling of GM and CSF (not required for WM) + Ylab{1} = Ylab{1} .* T3th(2) / cat_stat_kmeans( Ylab{1}(Ygm(:)>0.5), 1); + Ylab{3} = Ylab{3} .* T3th(1) / cat_stat_kmeans( Ylab{3}(Ycm(:)>0.5), 1); + + % RD202110: corrected cases with incorrect background region Ynb + if sum(Ynb(:)>0.5)>0 + Ynb = smooth3( Ycls{6})>128 & Yg<0.3; + end + if sum(Ynb(:)>0.5)>0 + Ylab{6} = Ylab{6} .* cat_stat_kmeans( Ysrc(Ynb(:)>0.5) , 1) / max(eps,cat_stat_kmeans( Ylab{6}(Ynb(:)>0.5) , 1)); + else + Ylab{6} = min(Ysrc(:)); + end + + %% restore original resolution if exist('resT0','var') Ycls = Yclso2; clear Yclso2; Ysrc = Ysrco2; clear Ysrco2; for i=[1 2 3 6], Ylab{i} = cat_vol_resize (Ylab{i} , 'dereduceV' , resT0 ); end - % for i=1:6, Ycls{i} = cat_vol_resize( Ycls{i} , 'dereduceV' , resT0 ); end - Yb = cat_vol_resize( single(Yb) , 'dereduceV' , resT0 )>0.5; + Yb2 = cat_vol_resize( single(Yb2) , 'dereduceV' , resT0 )>0.5; end - - %% local intensity modification of the original image % -------------------------------------------------------------------- - stime = cat_io_cmd(' Intensity transformation','g5','',verb,stime); dispc=dispc+1; + cat_io_cmd(' Intensity transformation','g5','',verb,stime); Yml = zeros(size(Ysrc)); - Yml = Yml + ( (Ysrc>=Ylab{2} ) .* (3 + (Ysrc-Ylab{2}) ./ max(eps,Ylab{2}-Ylab{3})) ); - Yml = Yml + ( (Ysrc>=Ylab{1} & Ysrc=Ylab{3} & Ysrc=Ylab{2} ) .* (3 + (Ysrc - Ylab{2}) ./ max(eps,Ylab{2} - Ylab{3})) ); + Yml = Yml + ( (Ysrc>=Ylab{1} & Ysrc=Ylab{3} & Ysrc10)=10; - %% global - if Tth.T3th(Tth.T3thx==1)==Tth.T3thx(Tth.T3thx==1) %invers + + %% update of the global intensity normalized map + if Tth.T3th(Tth.T3thx==1) == Tth.T3thx(Tth.T3thx==1) % invers intensities (T2/PD) Ymg = max(eps,(Ysrc + srcmin)./Ylab{2}); else - Ymg = (Ysrc + srcmin)./max(eps,(Ylab{2} + srcmin)); - Ymg = Ymg * Tth.T3th(Tth.T3thx==3)/(Tth.T3thx(Tth.T3thx==3)/3); + Ymg = Ysrc ./ max(eps,Ylab{2}); + Ymg = Ymg * Tthc.T3th(Tthc.T3thx==3)/(Tthc.T3thx(Tthc.T3thx==3)/3) + srcmin; % RD202004: corrected srcmin-correction end - clear Ylab Ysrc + if ~debug, clear Ylab Ysrc; end Ymg = cat_main_gintnorm(Ymg,Tth); - %% - if debug>1 - try %#ok % windows requires this... i don't know why - tpmci=tpmci+1; tmpmat = fullfile(pth,reportfolder,sprintf('%s_%s%02d%s.mat',nam,'LAS',tpmci,'postpeaks')); - save(tmpmat); - end - end + %% fill up CSF in the case of a skull stripped image if max(res.mn(res.lkp==5 & res.mg'>0.1)) < mean(res.mn(res.lkp==3 & res.mg'>0.3)) - YM = cat_vol_morph(Yb,'d'); + YM = cat_vol_morph(Yb2,'d'); Ymls = smooth3(max(Yml,YM*0.5)); - Yml(YM & Yml<0.5)=Ymls(YM & Yml<0.5); + Yml(YM & Yml<0.5) = Ymls(YM & Yml<0.5); clear Ymls YM end clear res - %% class correction and second logical class map Ycls2 + + + %% interpolate maps if an resolution adaptation was applied if exist('resT0','var') Ywm = cat_vol_resize( single(Ywm) , 'dereduceV' , resT0 )>0.5; Ygm = cat_vol_resize( single(Ygm) , 'dereduceV' , resT0 )>0.5; - Yp0 = cat_vol_resize( Yp0 , 'dereduceV' , resT0 ); end + + + + %% class correction and definition of the second logical class map Ycls2 Ynwm = Ywm & ~Ygm & Yml/3>0.95 & Yml/3<1.3; Ynwm = Ynwm | (smooth3(Ywm)>0.6 & Yml/3>5/6); Ynwm(smooth3(Ynwm)<0.5)=0; - Yngm = Ygm & ~Ywm & Yml/3<0.95; Yngm(smooth3(Yngm)<0.5)=0; - Yncm = ~Ygm & ~Ywm & ((Yml/3)>1/6 | Ycls{3}>128) & (Yml/3)<0.5 & Yb; - if ~debug, clear Ywm Ygm; else Yclso = Ycls; end + Yngm = Ygm & ~Ywm & Yml/3<0.95; Yngm(smooth3(Yngm)<0.5)=0; + Yncm = ~Ygm & ~Ywm & ((Yml/3)>1/6 | Ycls{3}>128) & (Yml/3)<0.5 & Yb2; + if debug==0, clear Ywm Ygm; end + + % cleanup outer GM that was mislabeled as WM + Yp0 = single( Ycls{1} )/255*2 + single( Ycls{2} )/255*3 + single( Ycls{3} )/255; % recreate label map Ycls{2} = cat_vol_ctype(single(Ycls{2}) + (Ynwm & ~Yngm & Yp0>=1.5)*256 - (Yngm & ~Ynwm & Yp0>=2)*256,'uint8'); Ycls{1} = cat_vol_ctype(single(Ycls{1}) - (Ynwm & ~Yngm & Yp0>=1.5)*256 + (Yngm & ~Ynwm & Yp0>=2)*256,'uint8'); %Ycls{3} = cat_vol_ctype(single(Ycls{3}) - ((Ynwm | Yngm) & Yp0>=2)*256,'uint8'); - %Ycls{3} = cat_vol_ctype(single(Ycls{3}) + (Yb & Yml<1.1 & ~Ynwm & ~Yngm)*256,'uint8'); - Ycls{1} = cat_vol_ctype(single(Ycls{1}) - (Yb & Yml<1.1 & ~Ynwm & ~Yngm)*256,'uint8'); - Ycls{2} = cat_vol_ctype(single(Ycls{2}) - (Yb & Yml<1.1 & ~Ynwm & ~Yngm)*256,'uint8'); + %Ycls{3} = cat_vol_ctype(single(Ycls{3}) + (Yb2 & Yml<1.1 & ~Ynwm & ~Yngm)*256,'uint8'); + Ycls{1} = cat_vol_ctype(single(Ycls{1}) - (Yb2 & Yml<1.1 & ~Ynwm & ~Yngm)*256,'uint8'); + Ycls{2} = cat_vol_ctype(single(Ycls{2}) - (Yb2 & Yml<1.1 & ~Ynwm & ~Yngm)*256,'uint8'); Ycls2 = {Yngm,Ynwm,Yncm}; - if ~debug, clear Yngm Ynwm Yncm Yb; end + if debug==0, clear Yngm Ynwm Yncm Yb2; end + + - %% + %% final correction of Yml similar to Ymg Yml = Yml/3; - %cat_io_cmd('','','',verb,stime); - % if debug - % cat_io_cmd(' ','','',verb,stime); - % else - % cat_io_cmd(' ','','',verb,stime); - % cat_io_cmd('cleanup',dispc,'',verb); - % end end diff --git a/cat_main_LAS2.m b/cat_main_LAS2.m deleted file mode 100644 index 8918d472..00000000 --- a/cat_main_LAS2.m +++ /dev/null @@ -1,902 +0,0 @@ -function [Yml,Ymg,Ycls,Ycls2,T3th] = ... - cat_main_LAS2(Ysrc,Ycls,Ym,Yb0,Yy,T3th,res,vx_vol,extopts,Tth) -% This is an exclusive subfunction of cat_main. -% ______________________________________________________________________ -% -% Local Adaptive Segmentation (LAS): -% -% This version of the local adaptive intensity correction includes a -% bias correction that is based on a maximum filter for the WM, a mean -% filter of GM and a mean filter of the head tissue (if available). -% For each tissue a refined logical map is generated to filter the local -% tissue intensity and approximate the local tissue intensity in the whole -% volume. Based on these values an intensity transformation is used. -% Compared to the global correction this has to be done for each voxel. -% To save time only a rough linear transformation is used. Finally, a -% second NLM-filter is used. -% It is important to avoid high intensity blood vessels in the process -% because they will push down local WM and GM intensity. -% -% ______________________________________________________________________ -% -% [Yml,Ymg,Ycls,Ycls2,T3th] = ... -% cat_main_LAS2(Ysrc,Ycls,Ym,Yb0,Yy,T3th,res,vx_vol,extopts,Tth) -% -% Yml .. local intensity normalized image -% Ymg .. global intensity normalized image -% Ycls .. corrected SPM tissue class map -% Ycls2 .. ? -% T3th .. tissue thresholds of CSF, GM, and WM in Ysrc -% -% Ysrc .. (bias corrected) T1 image -% Ycls .. SPM tissue class map -% Ym .. intensity corrected T1 image (BG=0,CSF=1/3,GM=2/3,WM=1) -% Yb0 .. brain mask -% Yy .. deformation map -% T3th .. intensity thresholds from global intensity normalization -% res .. SPM segmentation structure -% vx_vol .. voxel dimensions -% extopts .. -% Tth .. -% ______________________________________________________________________ -% -% Christian Gaser, Robert Dahnke -% Structural Brain Mapping Group (https://neuro-jena.github.io) -% Departments of Neurology and Psychiatry -% Jena University Hospital -% ______________________________________________________________________ -% $Id$ - - - -% ______________________________________________________________________ -% -% internal maps: -% -% Yg .. gradient map - edges between tissues -% Ydiv .. divergence map - sulci/gyri pattern, and blood vessels -% Yp0 .. label map - tissue classes (BG=0,CSF=1,GM=2,WM=3) -% -% Ysw .. save WM tissue map -% Ybv .. blood vessel map -% Ycp .. CSF / background area for distances estimation -% Ycd .. CSF / background distance map -% Ywd .. WM distance map -% -% Ycm .. CSF -% Ygm .. GM -% Ywm .. WM -% Yvt .. WM next to the ventricle map -% ______________________________________________________________________ -% -% Development / TODO: -% * replace old calls of cat_vol_morph by distance based operations with -% resolution parameter -% ______________________________________________________________________ - - - - % The reduction to 1 mm is not really good for the Ycls map. If this - % is required to support faster and memory saving preprocessing also - % for ultra-high resolution data further test are necessary! - %def.uhrlim = 0.7; - %extopts = cat_io_checkinopt(extopts,def); - extopts.uhrlim = 0.2; % no reduction for >0.4 mm - - - % set this variable to 1 for simpler debugging without reduceBrain - % function (that normally save half of processing time) - verb = extopts.verb - 1; - vxv = 1 / mean(vx_vol); % normalization of voxel size (mostly for old calls of cat_vol_morph) - dsize = size(Ysrc); - NS = @(Ys,s) Ys==s | Ys==s+1; % function to ignore brain hemisphere coding - LASstr = max(eps,min(1,extopts.LASstr)); % LAS strength (for GM/WM threshold) - manual correction based on R1109 (2017/02) - LAB = extopts.LAB; % atlas labels - LABl1 = 1; % use atlas map - cleanupstr = min(1,max(0,extopts.gcutstr)); % required to avoid critical regions (only used in case of atlas maps) - cleanupdist = min(3,max(1,1 + 2*cleanupstr)); - - - % set debug = 1 and do not clear temporary variables if there is a breakpoint in this file - dbs = dbstatus; debug = 0; - for dbsi = 1:numel(dbs), if strcmp(dbs(dbsi).name,'cat_main_LAS'); debug = 1; break; end; end - - - - -%% --------------------------------------------------------------------- -% First, we have to optimize the segments using further information that -% SPM do not use, such as the gradient, divergence and distance maps. -% (1) The gradient map (average of the first derivate of the T1 map) is -% an edge map and independent of the image intensity. It helps to -% avoid PVE regions and meninges. -% (2) The divergence (second derivate of the T1 map) help to identify -% sulcal and gyral pattern and therefore to find WM and CSF regions -% for further corrections and to avoid meninges and blood vessels. -% (3) Distance maps can help to describe the cortex by its expected -% thickness or to separate structures close to the skull or deep -% in the brain. -% -% Furthermore, special assumption can be used. -% The first is that the maximum property of the WM in T1 data that -% allows using of a maxim filter for the GM/WM region. -% The second is the relative stable estimation of CSF/BG that allows to -% estimate a distance map. Because, most regions have a thin layer of -% GM around the WM, we can avoid overestimation of the WM by the other -% maps (especially the divergence). -% --------------------------------------------------------------------- - - fprintf('\n'); - stime = cat_io_cmd(' Prepare maps','g5','',verb); - - - % general resolution limitation - % ------------------------------------------------------------------- - % Using a lot of maps make this function memory intensive that is - % critical for high resolution data. Furthermore, it can be expected - % that the full resolution is not required here. However, reducing - % the tissue class map can lead to other problems. Hence, the uhrlim - % parameter is maybe inactive (see above). - % ------------------------------------------------------------------- - if any( vx_vol < extopts.uhrlim/2 ) - % store old data that is needed in full resolution - Yclso2 = Ycls; - Ysrco2 = Ysrc; - - % reduce all input maps (Ysrc, Ycls, Yy, Ym, Yb0) - [Ysrc,resT0] = cat_vol_resize( Ysrc , 'reduceV' , vx_vol , extopts.uhrlim , 64 ); - for ci = 1:numel(Ycls) - Ycls{ci} = cat_vol_resize( Ycls{ci}, 'reduceV' , vx_vol , extopts.uhrlim , 64 ); - end - Yy2 = ones([size(Ysrc),4],'single'); - for ci = 1:numel(Yy) - Yy2(:,:,:,ci) = cat_vol_resize( Yy(:,:,:,ci) , 'reduceV' , vx_vol , extopts.uhrlim , 64 ); - end - Yy = Yy2; clear Yy2; - Ym = cat_vol_resize( Ym , 'reduceV' , vx_vol , extopts.uhrlim , 64 ); - Yb0 = cat_vol_resize( single(Yb0) , 'reduceV' , vx_vol , extopts.uhrlim , 64 )>0.5; - - vx_vol = resT0.vx_volr; - vxv = 1/ mean(resT0.vx_volr); - dsize = size(Ysrc); - end - - - % brain bounding box to save processing time - Yclso = Ycls; - [Ym,BB] = cat_vol_resize( Ym , 'reduceBrain' , vx_vol , round(10/mean(vx_vol)) , Yb0 ); - Yb = cat_vol_resize( Yb0 , 'reduceBrain' , vx_vol , round(10/mean(vx_vol)) , Yb0 ); clear Yb0; - - try - for i = 1:numel(Ycls), Ycls{i} = cat_vol_resize(Ycls{i} , 'reduceBrain' , vx_vol , BB.BB); end - catch - error('LAS failed probably due to bad tissue contrast.\n'); - end - - % helping maps (Yg = mean gradient = edge) and divergence (CSF or WM skeleton) - Yg = cat_vol_grad( Ym , vx_vol ); - Ydiv = cat_vol_div( max(0.33,Ym) , vx_vol ); - Yp0 = single( Ycls{1} )/255*2 + single( Ycls{2} )/255*3 + single( Ycls{3} )/255; % create label map - Yb = smooth3( Yb | cat_vol_morph(Yb,'d',2*vxv) & Ym<0.8 & Yg<0.3 & Ym>0 )>0.5; % increase brain mask, for missing GM - - - - - %% adding of atlas and template information (for subcortical structures) - % ------------------------------------------------------------------- -%%% 20180801 - Why not use the partitioning before LAS to benefit by the - % better ventricle, WMH and cerebellum maps? - % - Because the partitioning benefit in subcortical regions - % from LAS. >> call it twice? - % ------------------------------------------------------------------- - if LABl1 - stime = cat_io_cmd(' Prepare partitions','g5','',verb,stime); - - % map atlas to RAW space - for i=1:5 - try - Vl1A = spm_vol(extopts.cat12atlas{1}); - break - catch - % avoid read error in parallel processing - pause(rand(1)) - end - end - Yl1 = cat_vol_ctype( cat_vol_sample(res.tpm(1),Vl1A,Yy,0) ); - - - % load WM of the Dartel/Shooting Template for WMHs - use uint8 to save memory - % Ywtpm .. WM template map - Vtemplate = spm_vol(extopts.templates{end}); - Ywtpm = cat_vol_ctype( cat_vol_sample(res.tpm(1),Vtemplate(2),Yy,1) * 255,'uint8'); - if debug==0, clear Yy; end - Ywtpm = single(reshape(Ywtpm,dsize)); - spm_smooth(Ywtpm,Ywtpm,2*vxv); - Ywtpm = cat_vol_ctype(Ywtpm); - - - % apply boundary box for brain mask - Yl1 = cat_vol_resize(Yl1 ,'reduceBrain',vx_vol,round(4/mean(vx_vol)),BB.BB); - Ywtpm = cat_vol_resize(Ywtpm,'reduceBrain',vx_vol,round(4/mean(vx_vol)),BB.BB); - % The correction of image resolution is not required because of the adaptation of the Yy? - - - % do not reduce LASstr - LASmod = min(2,max(1,mean((Ym( NS(Yl1,LAB.BG) & Yg<0.1 & Ydiv>-0.05 & Ycls{1}>4)) - 2/3) * 8)); - else - LASmod = 1; %#ok - end - - - % adaptation of the LASstr depending on average basal values - LASstr = min(1,max(0.01,LASstr * LASmod)); % adaptation by local BG variation - LASfs = 1 / max(0.01,LASstr); % smoothing filter strength - LASi = min(8,round(LASfs)); % smoothing iteration (limited) - - - - - %% helping segments - % ------------------------------------------------------------------- - % We will now try to identify the CSF and WM regions those finally - % describe the GM area. - % * brain mask (Ybb), blood vessel maps (Ybv), - % * distance maps from CSF (Ycd), WM (Ywd), and brain mask (Ybd) - % * tissue maps for CSF (Ycp >> Ycm), GM (Ygm), and WM (Ysw >> Ywm) - % * ventricle map (Yvt) - % ------------------------------------------------------------------- - stime = cat_io_cmd(sprintf(' Prepare segments (LASmod = %0.2f)',LASmod),'g5','',verb,stime); - - % Ybb .. brain mask - % Don't trust SPM too much by using Yp0 because it may miss some GM areas! - % 20180801 - The brain mask should be correct now. - Ybb = cat_vol_morph((Yb & Ym>1.5/3 & Ydiv<0.05) | Yp0>1.5,'lo',vxv); - if debug==0, clear Yb; end - - - % Ysw .. save WM - % Ybv .. possible blood vessels - Ysw1 = cat_vol_morph(Ycls{2}>128 & (min(1,Ym)-Ydiv)<1.5,'lc',vxv*2) & (Ym - Ydiv)>5/6; - Ybv = ((min(1,Ym) - Ydiv + Yg)>2.0 | (Ycls{5}>16 & Ym<0.6 & Ycls{1}<192)) & ... - ~cat_vol_morph(Ysw1,'d',1) & Ym>0.2; - - - % Ycp = mask for CSF/BG distance initialization - Ycp1 = Ycls{3}>240 & Ydiv>0 & Yp0<1.1 & Ym<0.5 ; % typical CSF - Ycp2 = Ycls{5}>8 & Ycls{2}<32 & Ym<0.6 & Ydiv>0; % venes - Ycp3 = (Ym-Ydiv/4<0.4) & Ycls{3}>4 & Ycls{3}>16; % sulcal CSF - Ycp4 = (single(Ycls{6}) + single(Ycls{4}))>192; % non-CSF .. class 5 with error for negative t1 values - Ycp5 = ~cat_vol_morph(Ybb,'lc',5); % add background - Ycp6 = Ym<0.3; % but do not trust the brain mask! - Ycp = Ycp1 | Ycp2 | Ycp3 | Ycp4 | Ycp5 | Ycp6; % combine maps - Ycp(smooth3(Ycp)>0.4) = 1; % remove some meninges - if debug==0; clear Ycp1 Ycp2 Ycp3 Ycp4 Ycp5 Ycp6; end - - - % Ywd .. WM distance map - % Ysk .. divergence skeleton to improve CSF map - % deactivated due to problems in non-cortical structures - Ywd = cat_vbdist(single(Yp0>2.5),Yp0>0.5,vx_vol); % WM distance for skeleton - %Ysk = cat_vol_div(min(5,Ywd),2); % divergence skeleton - %Ysk = (Ym + min(0,Ysk))<0.2; % binary divergence skeleton - %Ycp = Ycp | Ysk; if debug==0, clear Ysk; end % combine skeleton with CSF map - - % Ycd .. CSF/BG distance map - Ycd = cat_vbdist(single(Ycp),~Ycp,vx_vol); % real CSF distance - Ycd((Ym-Ydiv<2/3 | Ydiv>0.1) & Ycls{3}>4 & Ycls{3}>1) = ... % correction for sulci - min(Ycd((Ym-Ydiv<2/3 | Ydiv>0.1) & Ycls{3}>4 & Ycls{3}>1),1.5); % ... maybe another distance estimation? - - - % we need to remove strong edge regions, because here is no GM layer between CSF and WM ??? - % Yb .. brain mask - % Ybd .. skull distance -%%% 20180801 - What is the difference of Yb2 to Ybb and the (previously replaced) Yb? - Yb2 = cat_vol_morph(~Ycp | (Ycls{3}>128),'lc',1); - Ybd = cat_vbdist(single(~Yb2),Yb2,vx_vol); - - - % Yvt .. Ventricle map WITHOUT atlas data as large deep CSF areas - % need the ventricle to avoid PVE GM between WM and CSF. - % There is an update for atlas data. - Yvt = (Yg + abs(Ydiv))>0.4 & smooth3(single(Ycls{1})/255)<0.5 & Ybd>20 & ... - cat_vol_morph(Ycls{3}>8,'d',vxv) & cat_vol_morph(Ycls{2}>8,'d',vxv); - Yvt = smooth3(Yvt)>0.7; - Yvt = smooth3(Yvt)>0.2; - - - - - %% final tissue maps: Ycm = CSF, Ygm = GM, Ywm = WM - % ------------------------------------------------------------------- - - % CSF: - Ycm = Ycp & Yb2 & Ycls{3}>192 & ~Ybv & Ym<0.45 & Yg<0.25 & Ym>0 ; - %Ycm = Ycm | (Yb2 & (Ym - max(0,Ydiv))<0.5); % adding of divergence information was not save - if debug==0, clear Ycp; end - - - %% WM: - % Different definitions of the possible WM (Ysw,Ysw2,Ysw3,Ysw4) - % and a variable (Ygw) to avoid non WM areas. - % Ysw1 .. defined above - % Ysw2 .. addition WM map that further include CSF distance (Ycd) and - % divergence information (Ydiv) and use a flexible intensity - % Ysw3 .. similar to Ysw2 with a skull-near and CSF distance criteria - % to reconstruct WM gyri - % Ysw4 .. was remove long ago - % Ygw .. map to of regions that we want to avoid in all possible WM - % definitions Ysw* - Ysw2 = Yb2 & (Ycd - Ydiv)>2 & Ydiv<0 & Ym>(0.9 + LASstr * 0.05); % general WM - Ysw3 = Yb2 & (Ycd - Ydiv .* Ycd)>4 & Ydiv<-0.01 & Ym>0.5 & Ybd<20 & Ycd>2; % addition WM gyri close to the cortex - %Ysw4 = ( (Ycd - Ydiv*5)>3 & Yb2 & Ym>0.4 & Ybd<20 & Ycd>2.5) & ... % further WM (rmoved long ago) - % (Ydiv<-0.01 & (Yg + max(0,0.05 - Ycd/100))<0.1 ); - - Ygw = ~Ybv & Yb2 & Ybd>1 & (Ycd>1.0 | (Yvt & Yp0>2.9)) & ... - (Yg + Ydiv<(Ybd/50) | (Ydiv - Ym)<-1); - - Ywm = ( Ycls{2}>252 | Ysw1 | Ysw2 | Ysw3 ) & Ygw; - if debug==0, clear Ysw Ysw2 Ysw3 Ygw; end - - - %% GM: - % Different criteria (Ygm1,Ygm2,Ygm3) to define the possible GM area - % whereas Ygm4 describe general GM limitation. - % Ygm1 .. upper intensity limitation - % Ygm2 .. lower intensity limitation - % Ygm3 .. avoid GM next to hard boundaries in the middle of the brain - Ygm1 = ( Ym - Ydiv - max(0,2 - Ycd)/10 )<0.9; % upper limit - Ygm2 = Ycls{1}>4 | ( Ym>0.7 & Ycls{3}>64 ) | Ycd<(Ym + Ydiv)*3; % lower limit - Ygm3 = smooth3(Yg>(Ybd/800) & Ycls{2}<240 )>0.6; % avoid GM-WM PVE - Ygm = Ybb & ~Yvt & ~Ybv & ~Ywm & ~Ycm & Ycd>0.5 & Ygm1 & Ygm2 & Ygm3; - if debug==0, clear Ybv Yvt Ygm1 Ygm2 Ygm3; end - - % Ygm4 .. general GM limitations - Ygm4 = Ybb & ~Ycm & ~Ywm & Ym>1/3 & Ym<2.8/3 & ... - Yg<0.4 & (Ym - Ydiv)>1/3 & (Ym - Ydiv)<1; - Ygm4( smooth3(Ygm4)<0.5 ) = 0; % remove small dots - - % combine possible area with limitation map - Ygm = Ygm | Ygm4; - if debug==0, clear Ygm4; end - Ygm = Ygm | (Ym>1.5/3 & Ym<2.8/3 & ~Ywm & ~Ycm & Ybb); - Ygm(smooth3(Ygm)<0.25) = 0; - if debug==0, clear Ybv Ycp; end - - - - - - %% further maps that require the atlas - if LABl1 - % SPM GM segmentation can be affected by inhomogeneities and some GM - % areas are miss-classified as CSF/GM (Ycls{5}) when the affine - % registration was not optimal. But for some regions we can trust - % these information. - - - - %% regions for cleanup of meninges and blood vessels - % Ycbp .. close to the cerebellum - % Ycbn .. not to deep in the cerebellum - % Ylhp / Yrhp .. GM next to the left/right hemisphere - % Ybv2 .. close to the skull and between the hemispheres or between cerebrum and cerebellum - % Ybvv .. sinus - Ycbp = cat_vbdist(single( NS(Yl1,LAB.CB)),Yb2,vx_vol); % close to cerebellum - Ycbn = cat_vbdist(single(~NS(Yl1,LAB.CB)),Yb2,vx_vol); % close to cerebellum surface - Ylhp = cat_vbdist(single(mod(Yl1,2)==1 & Yb2 & Yl1>0),Yb2,vx_vol); % GM next to the left hemisphere - Yrhp = cat_vbdist(single(mod(Yl1,2)==0 & Yb2 & Yl1>0),Yb2,vx_vol); % GM next to the right hemisphere - Ybv2 = Ycls{5}>2 & Ym<0.7 & Ym>0.3 & Yb2 & (... - ((Ylhp+Ybd/2)0.5; - Ybvv = (Ym - max(0,6 - abs(Ycbp - 6))/50)<0.6 & Ym>0.4 & Yb2 & Ycbp<8 & Ycbp>1; % sinus - if debug==0, clear Ycbp; end - - - - %% refinements of subcortical regions - % Thalamus (TH): - % YTH .. enlarged atlas ROI - % Ytd .. CSF distance in the thalamus YTH - % Yxd .. distance to the basal-ganglia (BGL) in the thalamus YTH - THth = 0.8 - LASstr*0.6; %0.5; % lower more thalamus - YTH = NS(Yl1,LAB.TH) | (cat_vol_morph(NS(Yl1,LAB.TH),'d',3) & Ym>0.5 & Ycls{1}>128); - Ytd = cat_vbdist(single(Ym<0.45),YTH | NS(Yl1,LAB.BG),vx_vol); Ytd(Ytd>2^16)=0; % CSF distance in the TH - Yxd = cat_vbdist(single(NS(Yl1,LAB.BG)),YTH,vx_vol); Yxd(Yxd>2^16)=0; % BGL distance in the TH - if debug==0, clear YTH; end - - - % Yss .. (enlarged) subcortical structures - % Ynw .. no WM ?? - Yss = NS(Yl1,LAB.BG) | NS(Yl1,LAB.TH); - Yss = Yss | (cat_vol_morph(Yss,'d',vxv*2) & Ym>2.25/3 & Ym<2.75/3 & Ydiv>-0.01); % add higher tissue around mask - Yss = Yss | (cat_vol_morph(Yss,'d',vxv*3) & NS(Yl1,LAB.VT) & Yp0>1.5 & Yp0<2.3); % add lower tissue around mask - Yss = Yss & Yp0>1.5 & (Yp0<2.75 | (Ym<(2.5 + LASstr * 0.45)/3 & Ydiv>-0.05)); % limit the enlarged map by intensity - Yss = Yss | ((Yxd./max(eps,Ytd + Yxd))>THth/2 & ... % save TH by distances - for overcorrected images - (Yp0<2.75 | (Ym<(2.75 + LASstr * 0.20)/3 & Ydiv>-0.05))); - Yss = cat_vol_morph(Yss,'o'); - - Ynw = (Yxd./max(eps,Ytd + Yxd))>THth/2 | (NS(Yl1,LAB.BG) & Ydiv>-0.01); - if debug==0, clear Ytd Yxd ; end - - - % increase CSF ROI - % Yvt .. improve ventricle ROI with atlas information -%%%%% 20180801 - Why not use the ventricle maps from the partitioning? - % - Because the partitioning is after LAS - % - Why not use a fast partitioning before LAS? - % Why not use the code from the partitioning? - % Is a good ventricle segmentation here not important? - Yvt = cat_vol_morph( (NS(Yl1,LAB.VT) | cat_vol_morph(Ycm,'o',3) ) ... - & Ycm & ~NS(Yl1,LAB.BG) & ~NS(Yl1,LAB.TH) & Ybd>30,'d',vxv*3) & ~Yss; - - - - %% meningeal corrections - % Ycx .. cerebellar and other CSF - % Ycwm .. cerebellar WM - % Yccm .. cerebellar CSF - % Ybwm .. brain WM - % Ybcm .. brain CSF - Ycx = ( NS(Yl1,LAB.CB) & ( (Ym - Ydiv)<0.55 | Ycls{3}>128 ) ) | ... - ( ( (Ym - Ydiv)<0.45 & Ycls{3}>8 ) | Ycls{3}>240 ); - % in the cerebellum tissue can be differentiated by div etc. - Ycwm = NS(Yl1,LAB.CB) & (Ym - Ydiv*4)>5/6 & Ycd>3 & Yg>0.05; - Yccm = NS(Yl1,LAB.CB) & Ydiv>0.02 & Ym<1/2 & Yg>0.05; - Ybwm = (Ym - Ydiv*4)>0.9 & Ycd>3 & Yg>0.05; - Ybcm = Ydiv>0.04 & Ym<0.55 & Yg>0.05; - - - - %% correction 1 of tissue maps - % Ywmtpm .. refined WM template map (original Ywtpm) - Ywmtpm = ( single( Ywtpm )/255 .* Ym .* (1 - Yg - Ydiv) .* ... - cat_vol_morph( NS(Yl1,1) .* Ybd/5 ,'e',1) )>0.6; % no WM hyperintensities in GM! - if debug==0, clear Ywtpm; end - - - % refine the GM map - Ygm = Ygm | (Yss & ~Yvt & ~Ycx & ~Ybv2 & ~Ycwm & ~(Yccm | Ybcm)); - Ygm = Ygm & ~Ywmtpm & ~Ybvv; % no WMH area - - - % refine the WM map - Ywm = Ywm & ~Yss & ~Ybv2 & ~Ynw; - Ywm = Ywm | Ycwm | Ybwm; - Ywmtpm(smooth3(Ywmtpm & Ym<11/12)<0.5) = 0; % 20180801 - Why here? - Ywm = Ywm & ~Ywmtpm & ~Ybvv & ~Yss; % no WM area - - - Ycm = Ycm | ( (Ycx | Yccm | Ybcm) & Yg<0.2 & Ym>0 & Ydiv>-0.05 & Ym<0.3 & Yb2 ) | Ybvv; - if debug==0, clear Ycwm Yccm Ycd Ybvv Ycx; end - - - % Mapping of the brain stem to the WM (well there were some small GM - % structures, but the should not effect the local segmentation to much. - % Ybs .. brain stem mask - Ybs = Ym<1.2 & Ym>0.9 & Yp0>1.5 & ... - cat_vol_morph(NS(Yl1,LAB.BS) & Ym<1.2 & Ym>0.9 & Yp0>2.5,'c',2*vxv); - Ygm = (Ygm & ~Ybs & ~Ybv2 & ~Ywm) | Yss; - Ywm = Ywm | (Ybs & Ym<1.1 & Ym>0.9 & Yp0>1.5) ; - end - - - - - - - - - %% back to original resolution for full bias field estimation - % ------------------------------------------------------------------- - Ycls = Yclso; clear Yclso; - Ycm = cat_vol_resize(Ycm , 'dereduceBrain' , BB); - Ygm = cat_vol_resize(Ygm , 'dereduceBrain' , BB); - Ywm = cat_vol_resize(Ywm , 'dereduceBrain' , BB); - - Yvt = cat_vol_resize(Yvt , 'dereduceBrain' , BB); - Yb2 = cat_vol_resize(Yb2 , 'dereduceBrain' , BB); - Yss = cat_vol_resize(Yss , 'dereduceBrain' , BB); - Ybb = cat_vol_resize(Ybb , 'dereduceBrain' , BB); - Ybs = cat_vol_resize(Ybs , 'dereduceBrain' , BB); - Ybv2 = cat_vol_resize(Ybv2 , 'dereduceBrain' , BB); - - Ym = cat_vol_resize(Ym , 'dereduceBrain' , BB); - Yp0 = cat_vol_resize(Yp0 , 'dereduceBrain' , BB); - Yl1 = cat_vol_resize(Yl1 , 'dereduceBrain' , BB); - Ybd = cat_vol_resize(Ybd , 'dereduceBrain' , BB); - Yg = cat_vol_resize(Yg , 'dereduceBrain' , BB); - Ydiv = cat_vol_resize(Ydiv , 'dereduceBrain' , BB); - Ywd = cat_vol_resize(Ywd , 'dereduceBrain' , BB); - - - % correction for negative values - [Ysrcx,thx] = cat_stat_histth(Ysrc,99); clear Ysrcx; %#ok - srcmin = thx(1); - Ysrc = Ysrc - srcmin; - T3th = T3th - srcmin; - Tthc = Tth; Tthc.T3th = Tth.T3th - srcmin; - if exist('Ysrco2','var'),Ysrco2 = Ysrco2 - srcmin; end - - - - - - - -%% --------------------------------------------------------------------- -% Now, we can estimate the local peaks -% --------------------------------------------------------------------- -% Estimation of the local WM intensity with help of intensity corrected -% GM, CSF and head tissues to avoid overfitting (eg. BWP cerebellum). -% CSF is problematic in high contrast or skull-stripped images and -% should not be used here or in the GM peak estimation. -% --------------------------------------------------------------------- - mres = 1.1; - stime = cat_io_cmd(' Estimate local tissue thresholds (WM)','g5','',verb,stime); - Ysrcm = cat_vol_median3(Ysrc.*Ywm,Ywm,Ywm); - rf = [10^5 10^4]; - T3th3 = max(1,min(10^6,rf(2) / (round(T3th(3)*rf(1))/rf(1)))); - Ysrcm = round(Ysrcm*T3th3)/T3th3; - - - % Ygwg .. large homogen GM regions (e.g. subcortical structures or in the BWP cerebellum) - % Ygwc .. large homogen CSF regions - Ygwg = Ycls{1}>128 & Ym>(2/3 - 0.04) & Ym<(2/3 + 0.04) & Ygm .*Ydiv>0.01; - Ygwg = Ygwg | (Ycls{1}>128 & Yg<0.05 & abs(Ydiv)<0.05 & ~Ywm & Ym<3/4); - Ygwc = Ycls{3}>128 & Yg<0.05 & ~Ywm & ~Ygm & Ywd<3; - Ygwc(smooth3(Ygwc)<0.5) = 0; - if debug==0, clear Ywd; end - - - % also use normalized GM tissues - % Ygi .. tissue-corrected intensity map (normalized for WM) - Ygi = Ysrc .* Ygwg * T3th(3)/mean(Ysrc(Ygwg(:))); - if cat_stat_nanmean(Ym(Ygwc))>0.1 - % use normalized CSF tissue only in images with clear CSF intensity - % (to avoid errors in skull-stripped data) - Ygi = Ygi + Ysrc .* Ygwc * T3th(3)/mean(Ysrc(Ygwc(:))); - end - - - %% limit image resolution for fast processing - [Yi,resT2] = cat_vol_resize( Ysrcm ,'reduceV',vx_vol,mres,32,'max'); % maximum reduction for the WM area - Ygi = cat_vol_resize( Ygi ,'reduceV',vx_vol,mres,32,'meanm'); % mean for other (intensity normalized) tissues - for xi = 1:2*LASi, Ygi = cat_vol_localstat(Ygi,Ygi>0,2,1); end; % intensity smoothing in the GM area - Ygi(smooth3(Ygi>0)<0.3) = 0; % remove small dots (meninges) - Yi = cat_vol_localstat(Yi,Yi>0,1,3); % maximum filtering to stabilize small WM structures - Yi(Yi==0 & Ygi>0) = Ygi(Yi==0 & Ygi>0); % combine the WM and GM map - for xi = 1:2*LASi, Yi = cat_vol_localstat(Yi,Yi>0,2,1); end % intensity smoothing in both regions no maximum here! - - - % Add head tissue if it is available. - if cat_stat_nanmean(Ym(Ygwc))>0.1 && cat_stat_nanmean(Ysrc(Ycls{6}(:)>128))>T3th(1) - try % no skull-stripped - % definion of head tissue - Ygwh = Ycls{6}>128 & ~Ygwc & ~Ygwg & Yg<0.5 & abs(Ydiv)<0.1 & ... - Ysrc>min( res.mn(res.lkp==6) * 0.5 ) & Ysrc ( cat_stat_nanmedian( Ysrc(Ygwh(:)) ) + 2 * cat_stat_nanstd( Ysrc(Ygwh(:)) ) ); - Ygwh( Ygwhn ) = false; clear Ygwhn - %% go to low resolution - [Ygi,resTB] = cat_vol_resize(Ysrc .* Ygwh * T3th(3)/mean(Ysrc(Ygwh(:))),'reduceV',vx_vol,mres,32,'meanm'); - % additional approximation - Ygi = cat_vol_approx(Ygi,'nh',resTB.vx_volr,2); Ygi = cat_vol_smooth3X(Ygi,2 * LASfs); - Ygi = Ygwh .* max(eps,cat_vol_resize(Ygi,'dereduceV',resTB)); - Yi(Yi==0) = Ygi(Yi==0); - if debug==0; clear Ygwh; end - end - end - if debug==0, clear Ygwg Ygwc; end - - - % final approximation of the WM inhomogeneity - Yi = cat_vol_approx(Yi,'nh',resT2.vx_volr,2); - Yi = cat_vol_smooth3X(Yi,2 * LASfs); - Ylab{2} = max(eps,cat_vol_resize(Yi,'dereduceV',resT2)); -% Ylab{2} = Ylab{2} .* mean( [median(Ysrc(Ysw(:))./Ylab{2}(Ysw(:))),1] ); - if debug==0; clear Ysw Yi; end - - - - - - - %% update GM tissue map after the WM bias correction - % -------------------------------------------------------------------- - Ygm(Ysrc./Ylab{2}>(T3th(2) + 0.90 * diff(T3th(2:3)))/T3th(3)) = 0; - Ygm(Ysrc./Ylab{2}<(T3th(2) + 0.75 * diff(T3th(2:3)))/T3th(3) & ... - Ysrc./Ylab{2}<(T3th(2) - 0.75 * diff(T3th(2:3)))/T3th(3) & ... - Ydiv<0.3 & Ydiv>-0.3 & Ybb & ~Ywm & ~Yvt & ~Ybv2 & Ycls{1}>48) = 1; - Ywd2 = cat_vbdist(single(Ywm),Yb2); % update WM distance - Ygmn = (Ywd2 - Ym + Ydiv)>0.5 & ~Ybv2 & ... - (Ym + 0.5 - Ydiv - Yg - Ywd2/10)>1/3 & ... low intensity tissue - ~(Ym - min(0.2,Yg + Ywd2/10 - Ydiv)<1/4) & ... - Yg0.5 & ~Ywm; - Ygmn(smooth3(Ygmn)<0.5) = 0; - Ygm = Ygm | Ygmn; % correct gm (SPM based) - if debug==0; clear Ygm4; end - - - - - - - %% update maps -%%% ... more details - - % Update CSF maps. - % Ycm .. updated CSF map - % Ycp1 .. typical CSF - % Ycp2 .. venes - % Ycp3 .. sulcal CSF - % Ycp4 .. save non-CSF - % Ycp5 .. but do not trust the brain mask! - % Ycpn .. updated CSF 2 map - % Ycd2 .. updated CSF distance - Ycm = Ycm | (Yb2 & (Ysrc./Ylab{2})<((T3th(1)*0.5 + 0.5*T3th(2))/T3th(3))); - Ycm = ~Ygm & ~Ywm & ~Ybv2 & Yg<0.6 & Ycm; - Ycp1 = Ycls{2}<128 & Ydiv>0 & Yp0<2.1 & Ysrc./Ylab{2}32 & Ycls{2}<32 & Ysrc./Ylab{2}0; % venes - Ycp3 = (Ym - Ydiv<0.4) & Ycls{3}>4 & Ycls{3}>16 & Ysrc./Ylab{2}192; % save non-CSF - Ycp5 = Ysrc./Ylab{2}0.4 ) = 1; % remove some meninges - Ycd2 = cat_vbdist(single(Ycpn),~Ycpn,vx_vol); - - - %% update GM map - Ygm = Ygm & ~Ycm & ~Ywm & Ywd2<5; - if debug==0; clear Ywd2; end - Ygm = Ygm | (NS(Yl1,1) & Ybd<20 & (Ycd2 - Ydiv)<2 & ... - Ycls{1}>0 & ~Ycm & Ybb & Ym>0.6 & Yg1.5 | ... - (Ym>1.2/3 & Ym<3.1/3 & Yb2))>0.6,'lo',min(1,vxv)); - if debug==0, clear Yp0 Yvt; end - Ygm(~Ybb) = 0; - Ygm(smooth3(Ygm)<0.3) = 0; - Ygm(smooth3(Ygm)>0.4 & Ysrc./Ylab{2}>mean(T3th(1)/T3th(2)) & ... - Ysrc./Ylab{2}<(T3th(2)*0.2+0.8*T3th(3))) = 1; - if debug==0; clear Ybb Yl1; end - - - - - - %% GM - % .. more details - - % Yi .. main GM intensity map - stime = cat_io_cmd(' Estimate local tissue thresholds (GM)','g5','',verb,stime); - Yi = Ysrc./Ylab{2} .* Ygm; - Yi = round(Yi * rf(2))/rf(2); - Yi(Ybs) = Ysrc(Ybs)./Ylab{2}(Ybs) .* T3th(2)/T3th(3); - if debug==0; clear Ybs; end - Yi = cat_vol_median3(Yi,Yi>0.5,Yi>0.5); % reduce artifacts - - % add CSF and CSF/GM areas (venes) to avoid overfitting - % Ycmx .. venes / sinus - % Tcmx .. intensity correction - Ycmx = smooth3( Ycm & Ysrc<(T3th(1)*0.8 + T3th(2)*0.2) )>0.9; - Tcmx = mean( Ysrc(Ycmx(:)) ./ Ylab{2}(Ycmx(:)) ) * T3th(3); - Yi(Ycmx) = Ysrc(Ycmx)./Ylab{2}(Ycmx) .* T3th(2)/Tcmx; - if debug==0; clear Ycm Ycmx; end - - %% reduce data to lower resolution - Ybgx = Ycls{6}>128; - Ybgxn = Ysrc < ( cat_stat_nanmedian( Ysrc(Ybgx(:)) ) - 2 * cat_stat_nanstd( Ysrc(Ybgx(:)) ) ) | ... - Ysrc > ( cat_stat_nanmedian( Ysrc(Ybgx(:)) ) + 2 * cat_stat_nanstd( Ysrc(Ybgx(:)) ) ); - Ybgx = Ybgx & ~Ybgxn; clear Ybgxn; - [Yi,resT2] = cat_vol_resize( Yi ,'reduceV',vx_vol,1,32,'meanm'); - Yii = cat_vol_resize( Ylab{2}/T3th(3) ,'reduceV',vx_vol,1,32,'meanm'); - Ybgx = cat_vol_resize( Ybgx ,'reduceV',vx_vol,1,32,'meanm'); - for xi = 1:2*LASi, Yi = cat_vol_localstat(Yi,Yi>0,3,1); end % local smoothing - Yi = cat_vol_approx(Yi,'nh',resT2.vx_volr,2); % approximation of the GM bias - Yi = min(Yi,Yii*(T3th(2) + 0.90*diff(T3th(2:3)))/T3th(3)); % limit upper values to avoid wholes in the WM - Yi(Ybgx) = Yii(Ybgx) * cat_stat_nanmean(Yi(~Ybgx(:))); % use WM bias field in the background - if debug==0; clear Yii Ybgx; end - Yi = cat_vol_smooth3X(Yi,LASfs); % final low res smoothing - Ylab{1} = cat_vol_resize(Yi,'dereduceV',resT2).*Ylab{2}; % back to original resolution - if debug==0; clear Yi; end - - - % update CSF map - Ycm = (single(Ycls{3})/255 - Yg*4 + abs(Ydiv)*2)>0.5 & ... - Ysrc<(Ylab{1}*mean(T3th([1,1:2]))/T3th(2)); - if debug==0; clear Ydiv; end - Ycm(smooth3(Ycm)<0.5)=0; - Ycm(Yb2 & cat_vol_morph(Ysrc128 & Yg<0.15; % & Ysrc0.5) * rf(2))/rf(2); - Yx = round(Ysrc./Ylab{2} .* Ynb * rf(2))/rf(2); - - % The correction should be very smooth (and fast) for CSF and we can use - % much lower resolution. - [Yc,resT2] = cat_vol_resize(Yc,'reduceV',vx_vol,8,16,'min'); % only pure CSF !!! - Yx = cat_vol_resize(Yx,'reduceV',vx_vol,8,16,'meanm'); - - % correct Nans and remove outlier - Yc(isnan(Yc)) = 0; Yc = cat_vol_median3(Yc,Yc~=0,true(size(Yc)),0.1); - Yx(isnan(Yx)) = 0; Yx = cat_vol_median3(Yx,Yx~=0,true(size(Yc)),0.1); - - % CSF intensity show be higher than background intensity - Yx(Yc>0)=0; Yc(Yx>0)=0; - if 0 % no so important ... keep it simple? - if cat_stat_nanmean(Yx(Yx~=0)) < cat_stat_nanmean(Yc(Yc~=0)) - meanYx = min(median(Yc(Yc(:)>0)),median(Yx(Yx(:)>0))); - meanYc = max(median(Yc(Yc(:)>0)),median(Yx(Yx(:)>0))); - else - meanYc = min(median(Yc(Yc(:)>0)),median(Yx(Yx(:)>0))); - meanYx = max(median(Yc(Yc(:)>0)),median(Yx(Yx(:)>0))); - end - else - meanYx = cat_stat_nanmean( [ 0 ; Yx(Yx(:)>0) ]); - meanYc = cat_stat_nanmedian( [ T3th(1) ; Yc(Yc(:)>0) ]); - end - - % correct for high intensity backgrounds (e.g. Gabi's R1) - %Yx = Yx .* (Yx > (meanYc + max(0.1,std(Yc(Yc(:)>0))))); - - %% avoid CSF in background and background in CSF - if ~res.ppe.affreg.highBG - bcd = max( abs( [ diff( [meanYx meanYc] ) min( diff( T3th/max(T3th(:)) ) ) ] )); - Yxc = Yx + (Yx==0 & Yc~=0) .* (Yc * meanYx/meanYc - 0.5 * bcd) + ... % avoid BG in the CSF (white dots in the blue CSF) - ( (Yx~=0) * 0.2 * bcd ); % this term is the lower HG intensity (some nice noise in the BG) - if 0% (meanYx/meanYc) > 0.8 % this looks nice but its no good and change thickness - Ycc = Yc + (Yc==0 & Yx~=0) .* (Yx * meanYc/meanYx + 0.5 * bcd); % here we avoid CSF in the BG (higher boudnary) ... blue CSF dots in the BG - else - Ycc = Yc; - end - %Ycc = Yc + min(max( meanYx + stdYbc , meanYc - stdYbc ), Yx .* meanYc/max(eps,meanYx)); - else - Yxc = Yx; - Ycc = Yc; - end - - %% guaranty small difference between backgound and CSF intensity - Yxa = cat_vol_approx(Yxc ,'nh',resT2.vx_volr,16); - Yca = cat_vol_approx(Ycc ,'nh',resT2.vx_volr,16); - if debug==0, clear Yc Yx; end - Yca = Yca * 0.7 + 0.3 * max( mean(Yca(:)) , T3th(1)/T3th(3) ); - % smoothing - Yxa = cat_vol_smooth3X( Yxa , LASfs * 2 ); - Yca = cat_vol_smooth3X( Yca , LASfs * 2 ); - %% back to original resolution, combination with WM bias, and final smoothing - Yca = cat_vol_resize(Yca,'dereduceV',resT2) .* Ylab{2}; - Yxa = cat_vol_resize(Yxa,'dereduceV',resT2) .* Ylab{2}; - Ylab{3} = cat_vol_smooth3X( Yca , LASfs * 2 ); - Ylab{6} = cat_vol_smooth3X( Yxa , LASfs * 2 ); - if debug==0, clear Yxa Yca; end - - if res.ppe.affreg.highBG - Ylab{6} = min(Ysrc(:)); - end - else - % simple - Ynb = smooth3( Ycls{6})>128 & Yg<0.3; - %Ynb = cat_vol_morph(Ynb,'e',4*vxv); - - Ylab{3} = ones(size(Ygm),'single') .* min( [ T3th(1) , cat_stat_nanmean(Ysrc(Ycm(:))) ] ); - Ylab{6} = ones(size(Ynb),'single') .* min( [ T3th(1) - min(diff(T3th)) , cat_stat_nanmean(Ysrc(Ycm(:))) ] ); - end - - % final scaling of GM and CSF (not required for WM) - Ylab{1} = Ylab{1} .* T3th(2) / cat_stat_kmeans( Ylab{1}(Ygm(:)>0.5), 1); - Ylab{3} = Ylab{3} .* T3th(1) / cat_stat_kmeans( Ylab{3}(Ycm(:)>0.5), 1); - - % RD202110: corrected cases with incorrect background region Ynb - if sum(Ynb(:)>0.5)>0 - Ynb = smooth3( Ycls{6})>128 & Yg<0.3; - end - if sum(Ynb(:)>0.5)>0 - Ylab{6} = Ylab{6} .* cat_stat_kmeans( Ysrc(Ynb(:)>0.5) , 1) / max(eps,cat_stat_kmeans( Ylab{6}(Ynb(:)>0.5) , 1)); - else - Ylab{6} = min(Ysrc(:)); - end - - %% restore original resolution - if exist('resT0','var') - Ycls = Yclso2; clear Yclso2; - Ysrc = Ysrco2; clear Ysrco2; - for i=[1 2 3 6], Ylab{i} = cat_vol_resize (Ylab{i} , 'dereduceV' , resT0 ); end - Yb2 = cat_vol_resize( single(Yb2) , 'dereduceV' , resT0 )>0.5; - end - - %% local intensity modification of the original image - % -------------------------------------------------------------------- - cat_io_cmd(' Intensity transformation','g5','',verb,stime); - - Yml = zeros(size(Ysrc)); - Yml = Yml + ( (Ysrc>=Ylab{2} ) .* (3 + (Ysrc - Ylab{2}) ./ max(eps,Ylab{2} - Ylab{3})) ); - Yml = Yml + ( (Ysrc>=Ylab{1} & Ysrc=Ylab{3} & Ysrc10)=10; - - - - %% update of the global intensity normalized map - if Tth.T3th(Tth.T3thx==1) == Tth.T3thx(Tth.T3thx==1) % invers intensities (T2/PD) - Ymg = max(eps,(Ysrc + srcmin)./Ylab{2}); - else - Ymg = Ysrc ./ max(eps,Ylab{2}); - Ymg = Ymg * Tthc.T3th(Tthc.T3thx==3)/(Tthc.T3thx(Tthc.T3thx==3)/3) + srcmin; % RD202004: corrected srcmin-correction - end - if ~debug, clear Ylab Ysrc; end - Ymg = cat_main_gintnorm(Ymg,Tth); - - - - %% fill up CSF in the case of a skull stripped image - if max(res.mn(res.lkp==5 & res.mg'>0.1)) < mean(res.mn(res.lkp==3 & res.mg'>0.3)) - YM = cat_vol_morph(Yb2,'d'); - Ymls = smooth3(max(Yml,YM*0.5)); - Yml(YM & Yml<0.5) = Ymls(YM & Yml<0.5); - clear Ymls YM - end - clear res - - - - %% interpolate maps if an resolution adaptation was applied - if exist('resT0','var') - Ywm = cat_vol_resize( single(Ywm) , 'dereduceV' , resT0 )>0.5; - Ygm = cat_vol_resize( single(Ygm) , 'dereduceV' , resT0 )>0.5; - end - - - - %% class correction and definition of the second logical class map Ycls2 - Ynwm = Ywm & ~Ygm & Yml/3>0.95 & Yml/3<1.3; - Ynwm = Ynwm | (smooth3(Ywm)>0.6 & Yml/3>5/6); Ynwm(smooth3(Ynwm)<0.5)=0; - Yngm = Ygm & ~Ywm & Yml/3<0.95; Yngm(smooth3(Yngm)<0.5)=0; - Yncm = ~Ygm & ~Ywm & ((Yml/3)>1/6 | Ycls{3}>128) & (Yml/3)<0.5 & Yb2; - if debug==0, clear Ywm Ygm; end - - % cleanup outer GM that was mislabeled as WM - Yp0 = single( Ycls{1} )/255*2 + single( Ycls{2} )/255*3 + single( Ycls{3} )/255; % recreate label map - Ycls{2} = cat_vol_ctype(single(Ycls{2}) + (Ynwm & ~Yngm & Yp0>=1.5)*256 - (Yngm & ~Ynwm & Yp0>=2)*256,'uint8'); - Ycls{1} = cat_vol_ctype(single(Ycls{1}) - (Ynwm & ~Yngm & Yp0>=1.5)*256 + (Yngm & ~Ynwm & Yp0>=2)*256,'uint8'); - %Ycls{3} = cat_vol_ctype(single(Ycls{3}) - ((Ynwm | Yngm) & Yp0>=2)*256,'uint8'); - %Ycls{3} = cat_vol_ctype(single(Ycls{3}) + (Yb2 & Yml<1.1 & ~Ynwm & ~Yngm)*256,'uint8'); - Ycls{1} = cat_vol_ctype(single(Ycls{1}) - (Yb2 & Yml<1.1 & ~Ynwm & ~Yngm)*256,'uint8'); - Ycls{2} = cat_vol_ctype(single(Ycls{2}) - (Yb2 & Yml<1.1 & ~Ynwm & ~Yngm)*256,'uint8'); - Ycls2 = {Yngm,Ynwm,Yncm}; - if debug==0, clear Yngm Ynwm Yncm Yb2; end - - - - %% final correction of Yml similar to Ymg - Yml = Yml/3; - -end - diff --git a/cat_main_LASs.m b/cat_main_LASs.m index 91b754bc..cf6ddb5c 100644 --- a/cat_main_LASs.m +++ b/cat_main_LASs.m @@ -279,6 +279,7 @@ else wmtht = 1; end + %% tissue thresholds % you need to know which are the extrem classes to define there PVE [gmm,gms,gmn] = cat_stat_kmeans(Ysrc(cat_vol_morph(Ygm & Ycls{1}>128,'e')),3); % we need 3 GM types to include GM/WM and GM/CSF PVE voxels @@ -288,56 +289,35 @@ TSmx = (TSth == max(TSth)) - (TSth == min(TSth)); TSth = TSth + TSmx .* [mean(gms(wmn>0.5)) mean(wms(wmn>0.5)) mean(cms(cmn>0.5))]; - %% - if 1 - Ysrcm = cat_vol_localstat(Ysrc .* Ywm,Ywm,1,wmtht); - Ysrcs = cat_vol_localstat(Ysrcm,Ysrcm>0,2,4); - Ysrcm(Ysrcs>mean(Ysrcs(Ysrcs(:)>0) + std(Ysrcs(Ysrcs(:)>0) ))) = 0; - Ysrcm = cat_vol_noPVE(Ysrcm,res.isMP2RAGE,vx_vol,2); - Ysrca = Ysrcm; - Ysrcmw = cat_vol_approx(Ysrcm,'nh',vx_vol,4); Ysrcmw = cat_vol_smooth3X(Ysrcmw,2); - - Ysrcm = cat_vol_localstat(Ysrc ./ TSth(1) .* TSth(2) .* (Ycls{1}>4),(Ygm & Ycls{1}>4),1,1);%Ygm & - Ysrcm = cat_vol_noPVE(Ysrcm,res.isMP2RAGE,vx_vol,2); - Ysrcs = cat_vol_localstat(Ysrcm,Ysrcm>0,4,4); - Ysrcm(Ysrcs>mean(Ysrcs(Ysrcs(:)>0)+std(Ysrcs(Ysrcs(:)>0) ))) = 0; - Ysrca = Ysrca + Ysrcm; - Ysrcmg = cat_vol_approx(Ysrcm,'nh',vx_vol,4); Ysrcmg = cat_vol_smooth3X(Ysrcmg,LASfs); - - Ysrcm = cat_vol_localstat(Ysrc ./ TSth(3) .* TSth(2) .* (Ycm & Ycls{3}>4),(Ycm & Ycls{3}>4),1,1); - Ysrcm = cat_vol_noPVE(Ysrcm,res.isMP2RAGE,vx_vol,2); - Ysrcs = cat_vol_localstat(Ysrcm,Ysrcm>0,2,4); - Ysrcm(Ysrcs>mean(Ysrcs(Ysrcs(:)>0)+std(Ysrcs(Ysrcs(:)>0) ))) = 0; - Ysrca = Ysrca + Ysrcm; - Ysrcmc = cat_vol_approx(Ysrcm,'nh',vx_vol,4); Ysrcmc = cat_vol_smooth3X(Ysrcmc,LASfs); - - Ysrcma = cat_vol_approx(Ysrca,'nh',vx_vol,2); Ysrcma = cat_vol_smooth3X(Ysrcma,LASfs); - Ysrcs = cat_vol_localstat(Ysrcma,Ysrcma>0,2,4); - Ysrcma(Ysrcs>mean(Ysrcs(Ysrcs(:)>0)+std(Ysrcs(Ysrcs(:)>0) ))) = 0; - - Yi = cat_vol_approx(Ysrcma,'nh',vx_vol,2); Yi = cat_vol_smooth3X(Yi,2); - Ylab{2} = Yi; - else - - Ysrcm = cat_vol_localstat(Ysrc .* Ywm,Ywm,1,wmtht) + ... - cat_vol_localstat(Ysrc ./ TSth(1) .* TSth(2) .* (Ygm & Ycls{1}>128),(Ygm & Ycls{1}>128),1,1) + ... - cat_vol_localstat(Ysrc ./ TSth(3) .* TSth(2) .* (Ycm & Ycls{3}>240),(Ycm & Ycls{3}>240),1,1); - Ysrcs = cat_vol_localstat(Ysrcm,Ysrcm>0,2,4); - Ysrcm(Ysrcs>mean(Ysrcs(Ysrcs(:)>0)+std(Ysrcs(Ysrcs(:)>0) ))) = 0; - Ysrcm2 = Ysrcm; - YM = cat_vol_morph(Ysrcm>0,'d') & Yp0>1 & Ysrcm==0; Ysrcm2(YM) = Ysrc(YM); - Ysrcm2 = cat_vol_localstat(Ysrcm2,Ysrcm2>0,2,wmtht); Ysrcm(YM) = Ysrcm2(YM); - Ysrcm = cat_vol_noPVE(Ysrcm,res.isMP2RAGE,vx_vol,2); - - % major correction outside the brain - [Yi ,resT2] = cat_vol_resize(Ysrcm,'reduceV',vx_vol,mres,32,'meanm'); % maximum reduction for the WM - if ~debug, clear Ysrcm Ydiv; end - Yi(smooth3(Yi>0)<0.5)=0; Yi(smooth3(Yi>0)<0.5)=0; - - % Yi(Yi==0 & Ygi>0)=Ygi(Yi==0 & Ygi>0); if ~debug, clear Ygi; end - Yi = cat_vol_approx(Yi,'nh',resT2.vx_volr,4); Yi = cat_vol_smooth3X(Yi,LASfs); - Ylab{2} = max(eps,cat_vol_resize(Yi,'dereduceV',resT2)); - end + + %% WM + Ysrcm = cat_vol_localstat(Ysrc .* Ywm,Ywm,1,wmtht); + Ysrcs = cat_vol_localstat(Ysrcm,Ysrcm>0,2,4); + Ysrcm(Ysrcs>mean(Ysrcs(Ysrcs(:)>0) + std(Ysrcs(Ysrcs(:)>0) ))) = 0; + Ysrcm = cat_vol_noPVE(Ysrcm,res.isMP2RAGE,vx_vol,2); + Ysrca = Ysrcm; + Ysrcmw = cat_vol_approx(Ysrcm,'nh',vx_vol,4); Ysrcmw = cat_vol_smooth3X(Ysrcmw,2); + + Ysrcm = cat_vol_localstat(Ysrc ./ TSth(1) .* TSth(2) .* (Ycls{1}>4),(Ygm & Ycls{1}>4),1,1);%Ygm & + Ysrcm = cat_vol_noPVE(Ysrcm,res.isMP2RAGE,vx_vol,2); + Ysrcs = cat_vol_localstat(Ysrcm,Ysrcm>0,4,4); + Ysrcm(Ysrcs>mean(Ysrcs(Ysrcs(:)>0)+std(Ysrcs(Ysrcs(:)>0) ))) = 0; + Ysrca = Ysrca + Ysrcm; + Ysrcmg = cat_vol_approx(Ysrcm,'nh',vx_vol,4); Ysrcmg = cat_vol_smooth3X(Ysrcmg,LASfs); + + Ysrcm = cat_vol_localstat(Ysrc ./ TSth(3) .* TSth(2) .* (Ycm & Ycls{3}>4),(Ycm & Ycls{3}>4),1,1); + Ysrcm = cat_vol_noPVE(Ysrcm,res.isMP2RAGE,vx_vol,2); + Ysrcs = cat_vol_localstat(Ysrcm,Ysrcm>0,2,4); + Ysrcm(Ysrcs>mean(Ysrcs(Ysrcs(:)>0)+std(Ysrcs(Ysrcs(:)>0) ))) = 0; + Ysrca = Ysrca + Ysrcm; + Ysrcmc = cat_vol_approx(Ysrcm,'nh',vx_vol,4); Ysrcmc = cat_vol_smooth3X(Ysrcmc,LASfs); + + Ysrcma = cat_vol_approx(Ysrca,'nh',vx_vol,2); Ysrcma = cat_vol_smooth3X(Ysrcma,LASfs); + Ysrcs = cat_vol_localstat(Ysrcma,Ysrcma>0,2,4); + Ysrcma(Ysrcs>mean(Ysrcs(Ysrcs(:)>0)+std(Ysrcs(Ysrcs(:)>0) ))) = 0; + + Yi = cat_vol_approx(Ysrcma,'nh',vx_vol,2); Yi = cat_vol_smooth3X(Yi,2); + Ylab{2} = Yi; if ~debug, clear Yi; end @@ -403,6 +383,7 @@ end clear Ynba Yca; + %% back to original resolution if any(resTb.vx_vol ~= resTb.vx_volr) Ysrc = Ysrco; clear Ysrco; @@ -411,7 +392,8 @@ Yp0 = cat_vol_resize(Yp0,'dereduceV',resTb); %#ok Ywm = cat_vol_resize(Ywm,'dereduceV',resTb); %#ok end - %% + + labmn = zeros(1,6); for i=[1 2 3 6], if ndims(Ylab{i})==3, labmn(i) = mean(Ylab{i}(Ycls{i}(:)>128)); else, labmn(i) = mean(Ylab{i}); end; end [labmn,lo] = sort(labmn,'descend'); %clear labmn; %#ok @@ -430,8 +412,6 @@ %% update of the global intensity normalized map % ######## RD202105: this is not optimal ######### Ymg = Ysrc ./ max(eps,Ylab{2}) * Tth.T3th(5); - %Ymg = Ysrc ./ max(eps,Ylab{2}); - %Ymg = Ymg * max( Tthc.T3th(Tthc.T3thx>0 & Tthc.T3thx<4) ) / max( max(Tthc.T3thx(Tthc.T3thx>0 & Tthc.T3thx<4))/3); % RD202004: corrected srcmin-correction Ymg = cat_main_gintnorm(Ymg,Tth); if ~debug, clear Ylab Ysrc; end diff --git a/cat_run_job.m b/cat_run_job.m index b108c237..48e7158a 100644 --- a/cat_run_job.m +++ b/cat_run_job.m @@ -43,7 +43,6 @@ function cat_run_job(job,tpm,subj) global cat_err_res; % for CAT error report cat_err_res.stime = clock; cat_err_res.cat_warnings = cat_io_addwarning('reset'); % reset warnings - stime = clock; stime0 = stime; % overall processing time @@ -85,18 +84,23 @@ function cat_run_job(job,tpm,subj) [pp,ff,ee,ex] = spm_fileparts(job.data{subj}); %#ok % sometimes we have to remove .nii from filename if files were zipped catlog = fullfile(pth,reportfolder,['catlog_' strrep(ff,'.nii','') '.txt']); - if exist(catlog,'file'), delete(catlog); end % write every time a new file, turn this of to have an additional log file + if exist(catlog,'file'), delete(catlog); end % write every time a new file, turn this off to have an additional log file % check if not another diary is already written that is not the default- or catlog-file. if ~strcmpi(spm_check_version,'octave') olddiary = spm_str_manip( get(0,'DiaryFile') , 't'); - usediary = ~isempty(strfind( olddiary , 'diary' )) | ~isempty(strfind( olddiary , 'catlog_' )); + usediary = ~isempty(strfind( olddiary , 'diary' )) || ~isempty(strfind( olddiary , 'catlog_' )); if usediary diary(catlog); diary on; else - cat_io_cprintf('warn',sprintf('External diary log is writen to "%s".\n',get(0,'DiaryFile'))); + cat_io_cprintf('warn',sprintf('External diary log is written to "%s".\n',get(0,'DiaryFile'))); end + else + % always use diary and don't check for old one for Octave + usediary = 1; + diary(catlog); + diary on; end % print current CAT release number and subject file @@ -261,7 +265,9 @@ function cat_run_job(job,tpm,subj) ofname = fullfile(pp,[ff ee]); nfname = fullfile(pp,mrifolder,['n' ff '.nii']); if strcmp(ee,'.nii') - copyfile(ofname,nfname,'f'); + if ~copyfile(ofname,nfname,'f') + spm('alert!',sprintf('ERROR: Check file permissions for folder %s.\n',fullfile(pp,mrifolder)),'',spm('CmdLine'),1); + end elseif strcmp(ee,'.img') V = spm_vol(job.channel(n).vols{subj}); Y = spm_read_vols(V); @@ -399,6 +405,7 @@ function cat_run_job(job,tpm,subj) best_vx = max( min(vx_vol) ,job.extopts.restypes.(restype)(1)); vx_voli = min(vx_vol ,best_vx ./ ((vx_vol > (best_vx + job.extopts.restypes.(restype)(2)))+eps)); case 'optimal' + %% aniso = @(vx_vol) (max(vx_vol) / min(vx_vol)^(1/3))^(1/3); % penetration factor volres = @(vx_vol) repmat( round( aniso(vx_vol) * prod(vx_vol)^(1/3) * 10)/10 , 1 , 3); % volume resolution optresi = @(vx_vol) min( job.extopts.restypes.(restype)(1) , max( median(vx_vol) , volres(vx_vol) ) ); % optimal resolution @@ -413,7 +420,7 @@ function cat_run_job(job,tpm,subj) if any( (vx_vol ~= vx_voli) ) stime = cat_io_cmd(sprintf('Internal resampling (%4.2fx%4.2fx%4.2fmm > %4.2fx%4.2fx%4.2fmm)',vx_vol,vx_voli)); - if 0 + if 1 imat = spm_imatrix(Vi.mat); Vi.dim = round(Vi.dim .* vx_vol./vx_voli); imat(7:9) = vx_voli .* sign(imat(7:9)); @@ -421,9 +428,13 @@ function cat_run_job(job,tpm,subj) Vn = spm_vol(job.channel(n).vols{subj}); cat_vol_imcalc(Vn,Vi,'i1',struct('interp',2,'verb',0,'mask',-1)); else + %% Small improvement for CAT12.9 that uses the cat_vol_resize function rather than the simple interpolation. + % However, postive effects only in case of strong reductions >2, ie. it is nearly useless. + jobr = struct(); jobr.data = {Vi.fname}; - jobr.interp = -2005; % spline with smoothing in case of downsampling - jobr.verb = 0; + jobr.interp = -3005; % spline with smoothing in case of downsampling; default without smoothing -5; + jobr.verb = debug; + jobr.lazy = 0; jobr.prefix = ''; jobr.restype.res = vx_voli; % use other resolution for test cat_vol_resize(jobr); @@ -438,24 +449,19 @@ function cat_run_job(job,tpm,subj) clear Vi Vn; - %% APP bias correction (APPs1 and APPs2) + %% Affine Preprocessing (APP) with SPM % ------------------------------------------------------------ % Bias correction is essential for stable affine registration % but also the following preprocessing. This approach uses the % SPM Unified segmentation for intial bias correction of the % input data with different FWHMs (low to high frequency) and % resolutions (low to high). - % As far as SPM finally also gives us a nice initial segmentation - % why not use it for a improved maximum based correction?! % ------------------------------------------------------------ - % 20181222 if ~strcmp(job.extopts.species,'human'), job.extopts.APP = 5; end - if job.extopts.APP==1 + if job.extopts.APP == 1 job.subj = subj; [Ym,Ybg,WMth] = cat_run_job_APP_SPMinit(job,tpm,ppe,n,... ofname,nfname,mrifolder,ppe.affreg.skullstripped); end - - end @@ -470,9 +476,22 @@ function cat_run_job(job,tpm,subj) obj.biasfwhm = job.opts.biasfwhm; obj.tpm = tpm; obj.reg = job.opts.warpreg; - obj.samp = job.opts.samp; - obj.tol = job.opts.tol; - obj.lkp = []; + obj.samp = job.opts.samp; % resolution of SPM preprocessing (def. 3, 1.5 as last highest TPM optimal level) + obj.tol = job.opts.tol; % stopping criteria for SPM iteration of outer/inner loops + obj.newtol = 1 + ( isfield(job,'useprior') && ~isempty(job.useprior) ); + % stopping criteria for outer (=tol) and inner loop: + % -1-old SPM (>9 iters, inner=tol), + % 0-old CAT more outer iterations (>19 iter, inner=tol), + % 1-new optimal/faster with additional AUC criteria to have SPM minimum iterations (>9 iters, inner=1e-2) + % 2-new accurate with additioal AUC criteria but CAT minimum iterations (>19 iter, outer=tol, inner=1e-4) - like 0 + obj.lkp = []; + if ~strcmp('human',job.extopts.species) + % RD202105: There are multiple problems in primates and increased + % accuracy is maybe better (eg. 0.5 - 0.66) + scannorm = 0.7; %prod(obj.image.dims .* vx_vol).^(1/3) / 20; % variance from typical field fo view to normalized parameter + obj.samp = obj.samp * scannorm; % normalize by voxel size + obj.fwhm = obj.fwhm * scannorm; + end if all(isfinite(cat(1,job.tissue.ngaus))) for k=1:numel(job.tissue) obj.lkp = [obj.lkp ones(1,job.tissue(k).ngaus)*k]; @@ -536,15 +555,17 @@ function cat_run_job(job,tpm,subj) VF.pinfo(1:2,:) = VF.pinfo(1:2,:)/spm_global(VF); VG.pinfo(1:2,:) = VG.pinfo(1:2,:)/spm_global(VG); - % APP step 1 rough bias correction + % APP step 1 rough bias correction and preparation of the affine + % registration % -------------------------------------------------------------- % Already for the rough initial affine registration a simple % bias corrected and intensity scaled image is required, because % large head intensities can disturb the whole process. % -------------------------------------------------------------- % ds('l2','',vx_vol,Ym, Yt + 2*Ybg,obj.image.private.dat(:,:,:)/WMth,Ym,60) - if job.extopts.APP == 1070 && ~ppe.affreg.highBG - stime = cat_io_cmd('APP: Rough bias correction'); + if job.extopts.APP == 1070 && ~ppe.affreg.highBG && ... + ( ~isfield(job,'useprior') || isempty(job.useprior) ) + stime = cat_io_cmd('Affine preprocessing (APP)'); try Ysrc = single(obj.image.private.dat(:,:,:)); [Ym,Yt,Ybg,WMth] = cat_run_job_APP_init1070(Ysrc,vx_vol,job.extopts.verb); %#ok @@ -569,8 +590,7 @@ function cat_run_job(job,tpm,subj) end end - if ~( job.extopts.setCOM && ~( isfield(job,'useprior') && ~isempty(job.useprior) ) && ... - ~ppe.affreg.highBG && strcmp(job.opts.affreg,'prior') ) + if ~( job.extopts.setCOM && ~( isfield(job,'useprior') && ~isempty(job.useprior) ) && ~ppe.affreg.highBG ) stime = cat_io_cmd('Affine registration','','',1,stime); end @@ -584,17 +604,17 @@ function cat_run_job(job,tpm,subj) Ym2 = (Ym2 - ths(1)) ./ diff(ths) .* (1 - Ybg); VF.dat(:,:,:) = cat_vol_ctype(Ym2 * 200,'uint8'); end - clear WI ths; + clear WI; % smoothing resa = obj.samp*2; % definine smoothing by sample size VF1 = spm_smoothto8bit(VF,resa); - VG1 = spm_smoothto8bit(VG,resa); %#ok + VG1 = spm_smoothto8bit(VG,resa); - elseif APP == 1 + elseif job.extopts.APP == 1 % APP by SPM - VF.dt = [spm_type('float32') spm_platform('bigend')]; - VF.dat(:,:,:) = single(Ym); + VF.dt = [spm_type('UINT8') spm_platform('bigend')]; + VF.dat(:,:,:) = cat_vol_ctype(Ym * 200,'uint8'); VF.pinfo = repmat([1;0],1,size(Ym,3)); % smoothing @@ -602,13 +622,12 @@ function cat_run_job(job,tpm,subj) VF1 = spm_smoothto8bit(VF,resa); VG1 = spm_smoothto8bit(VG,resa); - elseif job.extopts.setCOM && ~( isfield(job,'useprior') && ~isempty(job.useprior) ) && ... - ~ppe.affreg.highBG && strcmp(job.opts.affreg,'prior') + elseif job.extopts.setCOM && ~( isfield(job,'useprior') && ~isempty(job.useprior) ) && ~ppe.affreg.highBG % standard approach (no APP) with static resa value and no VG smoothing stime = cat_io_cmd('Coarse affine registration'); resa = 8; VF1 = spm_smoothto8bit(VF,resa); - VG1 = VG; %#ok + VG1 = VG; [Ym,Yt,Ybg,WMth] = APPmini(obj,VF,job.extopts.histth); else stime = cat_io_cmd('Skip initial affine registration due to high-intensity background','','',1); @@ -617,19 +636,24 @@ function cat_run_job(job,tpm,subj) end %% prepare affine parameter - aflags = struct('sep',obj.samp,'regtype','subj','WG',[],'WF',[],'globnorm',1); % job.opts.affreg + aflags = struct('sep',obj.samp,'regtype','subj','WG',[],'WF',[],'globnorm',1); aflags.sep = max(aflags.sep,max(sqrt(sum(VG(1).mat(1:3,1:3).^2)))); aflags.sep = max(aflags.sep,max(sqrt(sum(VF(1).mat(1:3,1:3).^2)))); % use affine transformation of given (average) data for longitudinal mode - if isfield(job,'useprior') && ~isempty(job.useprior) && strcmp(job.opts.affreg,'prior') + if isfield(job,'useprior') && ~isempty(job.useprior) + % even in the development pipeline the prior is a good start ! priorname = job.useprior{1}; [pp,ff,ee,ex] = spm_fileparts(priorname); %#ok catxml = fullfile(pp,reportfolder,['cat_' ff '.xml']); % check that file exists and get affine transformation if exist(catxml,'file') - cat_io_cmd(sprintf(' Use affine transformation from:\n"%s"',ff),'','',1,stime); + if strcmp(job.opts.affreg,'prior') + fprintf('\nUse affine transformation from:\n%s\n',priorname); + else + fprintf('\nInitialize with affine transformation from:\n%s\n',priorname); + end stime = cat_io_cmd(' ',' ','',job.extopts.verb); xml = cat_io_xml(catxml); % sometimes xml file does not contain affine transformation @@ -640,11 +664,11 @@ function cat_run_job(job,tpm,subj) else Affine = xml.SPMpreprocessing.Affine; affscale = 1; - useprior = 1; + useprior = 1 + ~strcmp(job.opts.affreg,'prior'); end else - cat_io_addwarning([mfilename ':UseAffRegPrior'],... - sprintf('Affine prior file "%s" not found. \\\\nEstimate individual affine transformation. ',catxml),2,[1 2],catxml); + cat_io_cprintf('warn',sprintf('WARNING: File "%s" not found. Use individual affine transformation\n',catxml)); + Affine = eye(4); useprior = 0; end clear catxml; @@ -654,6 +678,8 @@ function cat_run_job(job,tpm,subj) obj.lkp(obj.lkp == 6) = []; obj.lkp = [ obj.lkp 6*ones(1,8) ]; else + %% + Affine = eye(4); useprior = 0; end @@ -711,29 +737,32 @@ function cat_run_job(job,tpm,subj) stime = cat_io_cmd('Affine registration','','',1,stime); if job.extopts.APP > 0 VF.dt = [spm_type('UINT8') spm_platform('bigend')]; - VF.pinfo(1:2,:) = VF.pinfo(1:2,:)/spm_global(VF); VF.pinfo = repmat([1;0],1,size(Ym,3)); VF.dat(:,:,:) = cat_vol_ctype(Ym*200); end - VF1 = spm_smoothto8bit(VF,aflags.sep); %#ok - VG1 = spm_smoothto8bit(VG,aflags.sep); %#ok + VF1 = spm_smoothto8bit(VF,aflags.sep); + VG1 = spm_smoothto8bit(VG,aflags.sep); try spm_plot_convergence('Init','Affine registration','Mean squared difference','Iteration'); catch spm_chi2_plot('Init','Affine registration','Mean squared difference','Iteration'); end + warning off + if ~exist('affscale','var'), affscale = 1.0; end evalc('[Affine1,affscale1] = spm_affreg(VG1, VF1, aflags, Affine, affscale);'); + warning on if ~any(any(isnan(Affine1(1:3,:)))) && affscale1>0.5 && affscale1<3, Affine = Affine1; end end clear VG1 VF1 else - Affine = eye(4); + % no affine registration and preprocessing at all and just prepare the data VF = spm_vol(obj.image(1)); [Ym,Yt,Ybg,WMth] = APPmini(obj,VF,job.extopts.histth); %#ok - %[Ym,Yt,Ybg,WMth] = cat_run_job_APP_init1070(single(obj.image.private.dat(:,:,:)),vx_vol,job.extopts.verb); %#ok if ~debug, clear Yt; end + useprior = 0; + Affine = eye(4); end @@ -747,6 +776,8 @@ function cat_run_job(job,tpm,subj) % these regions later. % We further discussed to use a separate mask images but finally decided % to keep this as simple as possible using no additional options! + % Moreover, we have to test here anyway to create warnings in case + % of inoptimal settings (e.g. no SLC but possible large lesions). obj.image0 = spm_vol(job.channel(1).vols0{subj}); Ysrc0 = spm_read_vols(obj.image0); Ylesion = single(Ysrc0==0 | isnan(Ysrc0) | isinf(Ysrc0)); clear Ysrc0; @@ -775,7 +806,6 @@ function cat_run_job(job,tpm,subj) end end - %% APP for spm_maff8 % optimize intensity range % we have to rewrite the image, because SPM reads it again @@ -854,9 +884,8 @@ function cat_run_job(job,tpm,subj) obj.msk.dat = uint8( Ymsk ); obj.msk = spm_smoothto8bit(obj.msk,0.1); end - clear Ymsk; end - clear Ysrc; + clear Ysrc Ymsk; else % defintion of basic variables in case of no APP obj.image.dat(:,:,:) = single(obj.image.private.dat(:,:,:)); @@ -872,10 +901,10 @@ function cat_run_job(job,tpm,subj) %% Fine affine Registration with automatic selection in case of multiple TPMs. % This may not work for non human data (or very small brains). % This part should be an external (coop?) function? - if useprior + if useprior==1 stime = cat_io_cmd('SPM preprocessing 1 (estimate 1 - use prior):','','',1,stime); elseif job.extopts.setCOM == 10 % no maffreg - stime = cat_io_cmd('SPM preprocessing 1 (estimate 1 - skip TPM registration):','','',1,stime); + stime = cat_io_cmd('SPM preprocessing 1 (estimate 1 - use no TPM registration):','','',1,stime); else stime = cat_io_cmd('SPM preprocessing 1 (estimate 1 - TPM registration):','','',1,stime); end @@ -884,11 +913,7 @@ function cat_run_job(job,tpm,subj) wo = warning('QUERY','MATLAB:RandStream:ActivatingLegacyGenerators'); wo = strfind( wo.state , 'on'); if wo, warning('OFF','MATLAB:RandStream:ActivatingLegacyGenerators'); end - if numel(job.opts.tpm)>1 - %% merging of multiple TPMs - obj2 = obj; obj2.image.dat(:,:,:) = max(0.0,Ym); - [Affine,obj.tpm,res0] = cat_run_job_multiTPM(job,obj2,Affine,ppe.affreg.skullstripped,1); %#ok - elseif strcmp(job.extopts.species,'human') + if strcmp(job.extopts.species,'human') %% only one TPM (old approach); spm_plot_convergence('Init','affine registration to TPM','Mean squared difference','Iteration'); @@ -978,7 +1003,7 @@ function cat_run_job(job,tpm,subj) %fliptest = 1; % 1 - test x>1, 2 - test for shearing %[ppe.affreg.flipped, ppe.affreg.flippedval,stime] = cat_vol_testflipping(obj,Affine,fliptest,stime,0); - if ~ppe.affreg.skullstripped + if isfield(ppe.affreg,'skullstripped') && ~ppe.affreg.skullstripped %% affreg with brainmask if debug [AffineN,Ybi,Ymi,Ym0] = cat_run_job_APRGs(Ym,Ybg,VF,Pb,Pbt,Affine,vx_vol,obj,job); %#ok @@ -1039,11 +1064,12 @@ function cat_run_job(job,tpm,subj) if numel(job.extopts.affmod)>6, job.extopts.affmod = job.extopts.affmod(1:6); end % remove too many if numel(job.extopts.affmod)<3, job.extopts.affmod(end+1:3) = job.extopts.affmod(1); end % isotropic if numel(job.extopts.affmod)<6, job.extopts.affmod(end+1:6) = 0; end % add translation + fprintf('\n Modify affine regitration (S=[%+3d%+3d%+3d], T=[%+3d%+3d%+3d])',job.extopts.affmod); sf = (100 - job.extopts.affmod(1:3)) / 100; imat = spm_imatrix(Affine); COMc = [eye(4,3), [ 0; -24 / mean(imat(7:9)); -12 / mean(imat(7:9)); 1] ]; imat = spm_imatrix(Affine * COMc); - imat(1:3) = imat(1:3) + job.extopts.affmod(4:6); + imat(1:3) = imat(1:3) - job.extopts.affmod(4:6); imat(7:9) = imat(7:9) .* sf; AffineMod = spm_matrix(imat) / COMc; @@ -1061,9 +1087,9 @@ function cat_run_job(job,tpm,subj) % ds('l2','a',0.5,Ysrc/WMth,Yb,Ysrc/WMth,Yb,140); warning off try - % inital estimate - stime = cat_io_cmd('SPM preprocessing 1 (estimate 2 - Unified segmentation):','','',job.extopts.verb-1,stime); - obj.tol = job.opts.tol; + %% inital estimate + stime = cat_io_cmd('SPM preprocessing 1 (estimate 2):','','',job.extopts.verb-1,stime); + obj.tol = job.opts.tol; % reset within loop % RD202012: Missclassification of GM as CSF and BG as tissue: % We observed problems with high-quality data (e.g. AVGs) and @@ -1131,6 +1157,9 @@ function cat_run_job(job,tpm,subj) end end end + if ~exist('res','var') + cat_io_printf('SPM preprocessing with default settings failed. Run backup settings. \n'); + end warning on; if job.opts.redspmres @@ -1139,14 +1168,16 @@ function cat_run_job(job,tpm,subj) end catch + %% cat_io_addwarning([mfilename ':ignoreErrors'],'Run backup function (IN DEVELOPMENT).',1,[1 1]); if isfield(obj.image,'dat') tmp = obj.image.dat; else tmp = spm_read_vols(obj.image); - dt2 = cat_io_strrep(spm_type(obj.image.dt(1)),{'float32','float64'},{'single','double'}); - obj.image.dat = eval(sprintf('%s(tmp);',dt2)); + dt2 = obj.image.dt(1); + dts = cat_io_strrep(spm_type(dt2),{'float32','float64'},{'single','double'}); + obj.image.dat = eval(sprintf('%s(tmp);',dts)); obj.image.pinfo = repmat([1;0],1,size(tmp,3)); end if exist('Ybi','var') @@ -1160,18 +1191,55 @@ function cat_run_job(job,tpm,subj) end suc = 0; - while obj.tol<1 - obj.tol = obj.tol * 10; + % try higher accuracy + while obj.tol>10e-9 && suc == 0 + obj.tol = obj.tol / 10; try res = cat_spm_preproc8(obj); suc = 1; end end + if suc == 0 + % try lower accuracy + while obj.tol<1 && suc == 0 + obj.tol = obj.tol * 10; + try + res = cat_spm_preproc8(obj); + suc = 1; + end + end + end if any( (vx_vol ~= vx_voli) ) || ~strcmp(job.extopts.species,'human') [pp,ff,ee] = spm_fileparts(job.channel(1).vols{subj}); delete(fullfile(pp,[ff,ee])); end + + if suc==0 + %% + mati = spm_imatrix(V.mat); + + error('cat_run_job:spm_preproc8',sprintf([ + 'Error in spm_preproc8. Check image and orientation. \n'... + ' Volume size (x,y,z): %8.0f %8.0f %8.0f \n' ... + ' Origin (x,y,z): %8.1f %8.1f %8.1f \n' ... + ' Rotation (deg): %8.1f %8.1f %8.1f \n' ... + ' Resolution: %8.1f %8.1f %8.1f \n'],... + V.dim,[mati(1:3),mati(4:6),mati(7:9),])); + end + + %% set internal image + if ~exist('dt2','var') + %tmp = obj.image.dat; + dt2 = obj.image.dt(1); + dts = cat_io_strrep(spm_type(dt2),{'float32','float64'},{'single','double'}); + end + obj.image.dat = eval(sprintf('%s(tmp);',dts)); + obj.image.pinfo = repmat([1;0],1,size(tmp,3)); + obj.image.dt(1) = dt2; + res.image.dat = eval(sprintf('%s(tmp);',dts)); + res.image.pinfo = repmat([1;0],1,size(tmp,3)); + res.image.dt(1) = dt2; end if ppe.affreg.skullstripped || job.extopts.gcutstr<0 % here we have to add manually our no variance background value of 0 @@ -1182,7 +1250,8 @@ function cat_run_job(job,tpm,subj) res.wp(end+1) = numel(res.wp) * eps; res.lkp(end+1) = 4; end - + warning on + if job.extopts.expertgui>1 %% print the tissue peaks mnstr = sprintf('\n SPM-US: ll=%0.6f, Tissue-peaks: ',res.ll); @@ -1198,10 +1267,6 @@ function cat_run_job(job,tpm,subj) cat_io_cmd(' ',' '); end fprintf('%5.0fs\n',etime(clock,stime)); - clear Ysrc; - - - %% check contrast (and convergence) % RD202006: SPM peak averaging @@ -1240,7 +1305,7 @@ function cat_run_job(job,tpm,subj) end - % updated tpm information for skull-stripped data should be available for cat_main + %% updated tpm information for skull-stripped data should be available for cat_main if isfield(obj.tpm,'bg1') && exist('ppe','var') && ( ppe.affreg.skullstripped || job.extopts.gcutstr<0 ) fname = res.tpm(1).fname; res.tpm = obj.tpm; @@ -1249,12 +1314,14 @@ function cat_run_job(job,tpm,subj) spm_progress_bar('Clear'); cat_progress_bar('Clear'); - %% call main processing + % call main processing res.tpm = obj.tpm.V; - if exist('ppe'), res.ppe = ppe; end res.stime = stime0; res.catlog = catlog; - res.Affine0 = res.Affine; %AffineMod; + res.Affine0 = res.Affine; + res.image0 = spm_vol(job.channel(1).vols0{subj}); + if exist('ppe','var'), res.ppe = ppe; end + if isfield(job.extopts,'affmod') && any(job.extopts.affmod) res.AffineUnmod = AffineUnmod; res.AffineMod = AffineMod; @@ -1264,7 +1331,7 @@ function cat_run_job(job,tpm,subj) % SPM segmentation in regions outside the TPM volume. res.bge = Ybge; end - res.image0 = spm_vol(job.channel(1).vols0{subj}); + if exist('Ylesion','var'), res.Ylesion = Ylesion; else, res.Ylesion = false(size(res.image.dim)); end; clear Ylesion; if exist('redspmres','var'); res.redspmres = redspmres; res.image1 = image1; end job.subj = subj; diff --git a/cat_run_job1639.m b/cat_run_job1639.m index 0084281d..1b82443c 100644 --- a/cat_run_job1639.m +++ b/cat_run_job1639.m @@ -872,132 +872,124 @@ function cat_run_job1639(job,tpm,subj) stime = cat_io_cmd('SPM preprocessing 1 (estimate 1 - TPM registration):','','',1,stime); end if ~isempty(job.opts.affreg) && useprior~=1 && job.extopts.setCOM ~= 10 % setcom == 10 - never use ... && strcmp('human',job.extopts.species) - if numel(job.opts.tpm)>1 - %% merging of multiple TPMs - obj2 = obj; obj2.image.dat(:,:,:) = max(0.0,Ym); - [Affine,obj.tpm,res0] = cat_run_job_multiTPM(job,obj2,Affine,ppe.affreg.skullstripped,1); %#ok - Affine3 = Affine; - elseif 1 %if strcmp(job.extopts.species,'human') - %% only one TPM (old approach); - spm_plot_convergence('Init','Fine affine registration','Mean squared difference','Iteration'); - warning off - - try - Affine2 = spm_maff8(obj.image(1),obj.samp,(obj.fwhm+1)*16,obj.tpm,Affine ,job.opts.affreg,80); - catch - Affine2 = spm_maff8(obj.image(1),obj.samp,(obj.fwhm+1)*16,obj.tpm,Affine ,job.opts.affreg); - end - scl1 = abs(det(Affine1(1:3,1:3))); - scl2 = abs(det(Affine2(1:3,1:3))); + spm_plot_convergence('Init','Fine affine registration','Mean squared difference','Iteration'); + warning off + + try + Affine2 = spm_maff8(obj.image(1),obj.samp,(obj.fwhm+1)*16,obj.tpm,Affine ,job.opts.affreg,80); + catch + Affine2 = spm_maff8(obj.image(1),obj.samp,(obj.fwhm+1)*16,obj.tpm,Affine ,job.opts.affreg); + end + scl1 = abs(det(Affine1(1:3,1:3))); + scl2 = abs(det(Affine2(1:3,1:3))); - if any(any(isnan(Affine2(1:3,:)))) - try - Affine2 = spm_maff8(obj.image(1),obj.samp,(obj.fwhm+1)*4,obj.tpm,Affine ,job.opts.affreg,80); - catch - Affine2 = spm_maff8(obj.image(1),obj.samp,(obj.fwhm+1)*4,obj.tpm,Affine ,job.opts.affreg); - end - if any(any(isnan(Affine2(1:3,:)))) - Affine2 = Affine; - end - else - % check for > 10% larger scaling - if ~strcmp(job.opts.affreg,'prior') && scl1 > 1.1*scl2 && job.extopts.setCOM ~= 11 % setcom == 11 - use always - stime = cat_io_cmd(' Use initial fine affine registration.','warn','',1,stime); - %fprintf('\n First fine affine registration failed.\n Use affine registration from previous step. '); - Affine2 = Affine1; - scl2 = scl1; - end - end + if any(any(isnan(Affine2(1:3,:)))) try - Affine3 = spm_maff8(obj.image(1),obj.samp,obj.fwhm,obj.tpm,Affine2,job.opts.affreg,80); + Affine2 = spm_maff8(obj.image(1),obj.samp,(obj.fwhm+1)*4,obj.tpm,Affine ,job.opts.affreg,80); catch - Affine3 = spm_maff8(obj.image(1),obj.samp,obj.fwhm,obj.tpm,Affine2,job.opts.affreg); + Affine2 = spm_maff8(obj.image(1),obj.samp,(obj.fwhm+1)*4,obj.tpm,Affine ,job.opts.affreg); end - - if ~any(any(isnan(Affine3(1:3,:)))) - scl3 = abs(det(Affine3(1:3,1:3))); - % check for > 5% larger scaling - if ~strcmp(job.opts.affreg,'prior') && scl2 > 1.05*scl3 && job.extopts.setCOM ~= 11 % setcom == 11 - use always - stime = cat_io_cmd(' Use previous fine affine registration.','warn','',1,stime); - %fprintf('\n Final fine affine registration failed.\n Use fine affine registration from previous step. '); - Affine = Affine2; - else - Affine = Affine3; - end - else % Affine3 failed, use Affine2 - Affine = Affine2; + if any(any(isnan(Affine2(1:3,:)))) + Affine2 = Affine; end - warning on else - Affine2 = Affine1; - Affine3 = Affine1; + % check for > 10% larger scaling + if ~strcmp(job.opts.affreg,'prior') && scl1 > 1.1*scl2 && job.extopts.setCOM ~= 11 % setcom == 11 - use always + stime = cat_io_cmd(' Use initial fine affine registration.','warn','',1,stime); + %fprintf('\n First fine affine registration failed.\n Use affine registration from previous step. '); + Affine2 = Affine1; + scl2 = scl1; + end end - - %% test for flipping - %fliptest = 2; - %[ppe.affreg.flipped, ppe.affreg.flippedval,stime] = cat_vol_testflipping(obj,Affine,fliptest,stime); - - if 0 - %% visual control for development and debugging - VFa = VF; VFa.mat = AffineMod * VF.mat; %Fa.mat = res0(2).Affine * VF.mat; - if isfield(VFa,'dat'), VFa = rmfield(VFa,'dat'); end - [Vmsk,Yb] = cat_vol_imcalc([VFa,spm_vol(Pb)],Pbt,'i2',struct('interp',3,'verb',0,'mask',-1)); - %[Vmsk,Yb] = cat_vol_imcalc([VFa;obj.tpm.V(1:3)],Pbt,'i2 + i3 + i4',struct('interp',3,'verb',0)); - %[Vmsk,Yb] = cat_vol_imcalc([VFa;obj.tpm.V(5)],Pbt,'i2',struct('interp',3,'verb',0)); - ds('d2sm','',1,Ym,Ym*0.5 + 0.5*Ym.*(Yb>0.5),round(size(Yb,3)*0.6)) + try + Affine3 = spm_maff8(obj.image(1),obj.samp,obj.fwhm,obj.tpm,Affine2,job.opts.affreg,80); + catch + Affine3 = spm_maff8(obj.image(1),obj.samp,obj.fwhm,obj.tpm,Affine2,job.opts.affreg); end - - - if isfield(ppe.affreg,'skullstripped') && ~ppe.affreg.skullstripped - %% affreg with brainmask - if debug - [Affine,Ybi,Ymi,Ym0] = cat_run_job_APRGs(Ym,Ybg,VF,Pb,Pbt,Affine,vx_vol,obj,job); %#ok + + if ~any(any(isnan(Affine3(1:3,:)))) + scl3 = abs(det(Affine3(1:3,1:3))); + % check for > 5% larger scaling + if ~strcmp(job.opts.affreg,'prior') && scl2 > 1.05*scl3 && job.extopts.setCOM ~= 11 % setcom == 11 - use always + stime = cat_io_cmd(' Use previous fine affine registration.','warn','',1,stime); + %fprintf('\n Final fine affine registration failed.\n Use fine affine registration from previous step. '); + Affine = Affine2; else - [Affine,Ybi] = cat_run_job_APRGs(Ym,Ybg,VF,Pb,Pbt,Affine,vx_vol,obj,job); + Affine = Affine3; end + else % Affine3 failed, use Affine2 + Affine = Affine2; end - - if ppe.affreg.skullstripped || job.extopts.gcutstr<0 - %% update number of SPM gaussian classes - Ybg = 1 - spm_read_vols(obj.tpm.V(1)) - spm_read_vols(obj.tpm.V(2)) - spm_read_vols(obj.tpm.V(3)); - noCSF = job.extopts.gcutstr == -2; - if 1 - for k=1:3 - noCSF - obj.tpm.dat{k} = spm_read_vols(obj.tpm.V(k)); - obj.tpm.V(k).dt(1) = 64; - obj.tpm.V(k).dat = double(obj.tpm.dat{k}); - obj.tpm.V(k).pinfo = repmat([1;0],1,size(Ybg,3)); - end + warning on + else + Affine2 = Affine1; + Affine3 = Affine1; + end + + %% test for flipping + %fliptest = 2; + %[ppe.affreg.flipped, ppe.affreg.flippedval,stime] = cat_vol_testflipping(obj,Affine,fliptest,stime); + + if 0 + %% visual control for development and debugging + VFa = VF; VFa.mat = AffineMod * VF.mat; %Fa.mat = res0(2).Affine * VF.mat; + if isfield(VFa,'dat'), VFa = rmfield(VFa,'dat'); end + [Vmsk,Yb] = cat_vol_imcalc([VFa,spm_vol(Pb)],Pbt,'i2',struct('interp',3,'verb',0,'mask',-1)); + %[Vmsk,Yb] = cat_vol_imcalc([VFa;obj.tpm.V(1:3)],Pbt,'i2 + i3 + i4',struct('interp',3,'verb',0)); + %[Vmsk,Yb] = cat_vol_imcalc([VFa;obj.tpm.V(5)],Pbt,'i2',struct('interp',3,'verb',0)); + ds('d2sm','',1,Ym,Ym*0.5 + 0.5*Ym.*(Yb>0.5),round(size(Yb,3)*0.6)) + end + + + if isfield(ppe.affreg,'skullstripped') && ~ppe.affreg.skullstripped + %% affreg with brainmask + if debug + [Affine,Ybi,Ymi,Ym0] = cat_run_job_APRGs(Ym,Ybg,VF,Pb,Pbt,Affine,vx_vol,obj,job); %#ok + else + [Affine,Ybi] = cat_run_job_APRGs(Ym,Ybg,VF,Pb,Pbt,Affine,vx_vol,obj,job); + end + end + + if ppe.affreg.skullstripped || job.extopts.gcutstr<0 + %% update number of SPM gaussian classes + Ybg = 1 - spm_read_vols(obj.tpm.V(1)) - spm_read_vols(obj.tpm.V(2)) - spm_read_vols(obj.tpm.V(3)); + noCSF = job.extopts.gcutstr == -2; + if 1 + for k=1:3 - noCSF + obj.tpm.dat{k} = spm_read_vols(obj.tpm.V(k)); + obj.tpm.V(k).dt(1) = 64; + obj.tpm.V(k).dat = double(obj.tpm.dat{k}); + obj.tpm.V(k).pinfo = repmat([1;0],1,size(Ybg,3)); end + end - obj.tpm.V(4 - noCSF).dat = Ybg; - obj.tpm.dat{4 - noCSF} = Ybg; - obj.tpm.V(4 - noCSF).pinfo = repmat([1;0],1,size(Ybg,3)); - obj.tpm.V(4 - noCSF).dt(1) = 64; - obj.tpm.dat(5 - noCSF:6) = []; - obj.tpm.V(5 - noCSF:6) = []; - obj.tpm.bg1(4 - noCSF) = obj.tpm.bg1(6); - obj.tpm.bg2(4 - noCSF) = obj.tpm.bg1(6); - obj.tpm.bg1(5 - noCSF:6) = []; - obj.tpm.bg2(5 - noCSF:6) = []; - %obj.tpm.V = rmfield(obj.tpm.V,'private'); - - % tryed 3 peaks per class, but BG detection error require manual - % correction (set 0) that is simple with only one class - % RD202306: SPM is not considering things without variation and - % a zeroed background is simply not existing! - % Moreover it is possible just to ignore classes :D - % Hence, we may not need to redefine the TPM at all. - if noCSF - job.opts.ngaus = [([job.tissue(1:2).ngaus])',1]; % gaussian background - else - job.opts.ngaus = ([job.tissue(1:3).ngaus])'; % no gaussian background - end - obj.lkp = []; - for k=1:numel(job.opts.ngaus) - job.tissue(k).ngaus = job.opts.ngaus(k); - obj.lkp = [obj.lkp ones(1,job.tissue(k).ngaus)*k]; - end + obj.tpm.V(4 - noCSF).dat = Ybg; + obj.tpm.dat{4 - noCSF} = Ybg; + obj.tpm.V(4 - noCSF).pinfo = repmat([1;0],1,size(Ybg,3)); + obj.tpm.V(4 - noCSF).dt(1) = 64; + obj.tpm.dat(5 - noCSF:6) = []; + obj.tpm.V(5 - noCSF:6) = []; + obj.tpm.bg1(4 - noCSF) = obj.tpm.bg1(6); + obj.tpm.bg2(4 - noCSF) = obj.tpm.bg1(6); + obj.tpm.bg1(5 - noCSF:6) = []; + obj.tpm.bg2(5 - noCSF:6) = []; + %obj.tpm.V = rmfield(obj.tpm.V,'private'); + + % tryed 3 peaks per class, but BG detection error require manual + % correction (set 0) that is simple with only one class + % RD202306: SPM is not considering things without variation and + % a zeroed background is simply not existing! + % Moreover it is possible just to ignore classes :D + % Hence, we may not need to redefine the TPM at all. + if noCSF + job.opts.ngaus = [([job.tissue(1:2).ngaus])',1]; % gaussian background + else + job.opts.ngaus = ([job.tissue(1:3).ngaus])'; % no gaussian background + end + obj.lkp = []; + for k=1:numel(job.opts.ngaus) + job.tissue(k).ngaus = job.opts.ngaus(k); + obj.lkp = [obj.lkp ones(1,job.tissue(k).ngaus)*k]; end end diff --git a/cat_run_job_multiTPM.m b/cat_run_job_multiTPM.m deleted file mode 100644 index cd9e6a1d..00000000 --- a/cat_run_job_multiTPM.m +++ /dev/null @@ -1,313 +0,0 @@ -function [Affine,tpm,res] = cat_run_job_multiTPM(job,obj,Affine,skullstripped,msk,acc) -% ______________________________________________________________________ -% Tissue Probability Maps (TPM) play an important role for the affine -% registration and segmentation. Preprocessing of atypical subjects -% (e.g. very old/young) with typical SPM TPM often lead to inoptimal -% results or can fail completelly. Hence, it is important to use at -% least specific TPMs in children. -% -% This function focus on the determination of the best fitting TPMs -% to create a subject specific TPM. It allows to create linear sub -% levels for age specific TPMs for e.g. very young subjects with -% small head, healty adults, and older subjects with severe brain -% atrophy. -% -% [Affine,tpm,res] = ... -% cat_run_job_multiTPM(job,obj[,Affine,skullstripped,msk,acc]) -% -% job .. CAT job structure -% obj .. SPM preprocessing structure -% Affine .. affine in/output matrix -% skullstripped .. use brain tissues only -% msk .. mask regions in the unified segmentation -% acc .. accuracy (default=0, just for tests) -% 0 - faster but less exact -% 1 - slower but more exact -% -% This function is part of the CAT preprocessing (cat_main_job) and -% utilize the TPM based affine registration spm_maffreg of SPM and the -% unified segmentation. -% ______________________________________________________________________ -% -% Christian Gaser, Robert Dahnke -% Structural Brain Mapping Group (https://neuro-jena.github.io) -% Departments of Neurology and Psychiatry -% Jena University Hospital -% ______________________________________________________________________ -% $Id$ -% ______________________________________________________________________ - -%#ok<*WNOFF,*WNON> - - dbs = dbstatus; debug = 0; for dbsi=1:numel(dbs), if strcmp(dbs(dbsi).name,mfilename); debug = 1; break; end; end - - stime0 = clock; - - if ~exist('Affine','var'), Affine = eye(4); end - if ~exist('skullstripped','var'), skullstripped = 0; end - if ~exist('acc','var') % overwrite parameter for tests - obj.samp = job.opts.samp + 1; % higher res (default: 3 mm) - obj.tol = min(1e-2,job.opts.tol * 10); % faster pp (default: 1-e4;) - else - acc = min(1,max(0,acc)); - obj.samp = 4 - 2*acc; % higher res (default: 3 mm) - obj.tol = 10^(-2 - 3*acc); % faster pp (default: 1-e4;) - end - if ~exist('msk','var'), msk = 0; end - - - % mask probably masked/stripped voxels! - if isfield(obj,'image0'), obj = rmfield(obj,'image0'); end - obj.image = cat_vol_resize(obj.image,'interpv',1.5); - if msk - if isfield(obj,'msk') - obj.msk = cat_vol_resize(obj.msk,'interpv',1.5); - else - obj.msk = obj.image(1); - obj.msk.dt = [spm_type('uint8') spm_platform('bigend')]; - obj.msk.pinfo = repmat([255;0],1,obj.image(1).dim(3)); - if numel(msk)==1 - obj.msk.dat = uint8(spm_read_vols(obj.image)==0); - else - obj.msk.dat = uint8(msk>0); - end - if isfield(obj.msk,'private'), obj.msk = rmfield(obj.msk,'private'); end - if isfield(obj.image,'private'), obj.image = rmfield(obj.image,'private'); end - obj.msk = spm_smoothto8bit(obj.msk,0.1); - end - else - if isfield(obj,'msk'), obj = rmfield(obj,'msk'); end - end - - - if skullstripped>1 - Vi = obj.image; - if isfield(Vi,'dat'), Vi = rmfield(Vi,'dat'); end - [Vmsk,Yb] = cat_vol_imcalc([Vi;obj.tpm.V(1:3)],'blub', ... - 'i2 + i3 + i4',struct('interp',5,'verb',0)); clear Vmsk; %#ok - obj.image.dat = obj.image.dat .* (Yb>0.1); - if ~debug, clear Yb; end - end - - - - % In case of skull-stripped images it is better to use a more simple - % model with only 4 classes with CSF, GM, WM and background. 3 peaks - % per tissue class seamed to be more stable as only 1 peak! - if skullstripped - job.opts.ngaus = 3*ones(4,1); % this is more safe - obj.lkp = []; - for k=1:numel(job.opts.ngaus) - job.tissue(k).ngaus = job.opts.ngaus(k); - obj.lkp = [obj.lkp ones(1,job.tissue(k).ngaus)*k]; - end - else - % simplyfied model? - %obj.lkp = 1:6; - end - - - % prepare variables - mAffine = cell(1,numel(job.opts.tpm)); - weight = zeros(1,numel(job.opts.tpm)); - - - % estimate affine registration to the current TPM - %res = struct(); - for i=1:numel(job.opts.tpm) - if numel(job.opts.tpm)>1 - spm_plot_convergence('Init',... - sprintf('Fine affine registration (TPM %d)',i),... - 'Mean squared difference','Iteration'); - if i==1 - fprintf('\n'); - stime = cat_io_cmd(sprintf(' Evaluate TPM %d (%s)',i,... - spm_str_manip(job.opts.tpm{i},'ra30')),'g5'); - else - stime = cat_io_cmd(sprintf(' Evaluate TPM %d (%s)',i,... - spm_str_manip(job.opts.tpm{i},'ra30')),'g5','',job.extopts.verb-1,stime); - end - else - stime = clock; - spm_plot_convergence('Init','Fine affine registration',... - 'Mean squared difference','Iteration'); - end - - - % select and load TPM - [pp,ff,ee] = spm_fileparts(job.opts.tpm{i}); - job.opts.tpm{i} = fullfile(pp,[ff,ee]); - obj.tpm = spm_load_priors8(job.opts.tpm{i}); - if all(numel(obj.tpm.V)~=[4 6]) - error('cat_run_job_multiTPM:badtpm',[... - 'The TPM has to include 4 (GM,WM,CSF,BG) or 6 (GM,WM,CSF,skull,head,BG)\n' ... - 'tissue classes! Your TPM has %d! Is this the right file?\n' ... - ' %s'],numel(obj.tpm.V),job.opts.tpm{i}); - end - - % update number of SPM gaussian classes - if skullstripped - Yb2 = 1 - spm_read_vols(obj.tpm.V(1)) - spm_read_vols(obj.tpm.V(2)) - spm_read_vols(obj.tpm.V(3)); - if 1 - for k=1:3 - obj.tpm.dat{k} = spm_read_vols(obj.tpm.V(k)); - obj.tpm.V(k).dt(1) = 64; - obj.tpm.V(k).dat = double(obj.tpm.dat{k}); - obj.tpm.V(k).pinfo = repmat([1;0],1,size(Yb2,3)); - end - end - - obj.tpm.V(4).dat = Yb2; - obj.tpm.dat{4} = Yb2; - obj.tpm.V(4).pinfo = repmat([1;0],1,size(Yb2,3)); - obj.tpm.V(4).dt(1) = 64; - obj.tpm.dat(5:6) = []; - obj.tpm.V(5:6) = []; - obj.tpm.bg1(4) = obj.tpm.bg1(6); - obj.tpm.bg2(4) = obj.tpm.bg1(6); - obj.tpm.bg1(5:6) = []; - obj.tpm.bg2(5:6) = []; - obj.tpm.V = rmfield(obj.tpm.V,'private'); - if ~debug, clear Yb2; end - end - - % affine registration to TPM with (rought) and without smoothing (fine) - try - warning off - mAffine{i} = spm_maff8(obj.image(1),... - obj.samp*2,(obj.fwhm+1)*10,obj.tpm,Affine,job.opts.affreg,160*(1+acc)); - mAffine{i} = spm_maff8(obj.image(1),... - obj.samp*2,(obj.fwhm+1)*5,obj.tpm,mAffine{i},job.opts.affreg,80*(1+acc)); - mAffine{i} = spm_maff8(obj.image(1),... - obj.samp,obj.fwhm,obj.tpm,mAffine{i},job.opts.affreg,40*(1+acc)); - mAffine{i} = spm_maff8(obj.image(1),... - obj.samp/2,obj.fwhm,obj.tpm,mAffine{i},job.opts.affreg,20*(1+acc)); - warning on - catch - mAffine{i} = Affine; - end - % further checks and corrections?? - - - % call Unified Segmentation - if ~any(isnan(mAffine{i}(:))) - obj.Affine = mAffine{i}; - try - res(i) = cat_spm_preproc8(obj); %#ok<*AGROW> - catch - res(i).wp = nan(1,6); - res(i).lkp = obj.lkp; - res(i).mn = nan(1,numel(obj.lkp)); - res(i).mg = nan(numel(obj.lkp),1); - res(i).vr = nan(1,1,numel(obj.lkp)); - end - else - res(i).wp = nan(1,6); - res(i).lkp = obj.lkp; - res(i).mn = nan(1,numel(obj.lkp)); - res(i).mg = nan(numel(obj.lkp),1); - res(i).vr = nan(1,1,numel(obj.lkp)); - end - - % estimate weighting - % The variable res(i).wp gives good information how good segmentation - % fits to the TPM. The first 3 classes should be relative similar and - % have high values (>0.1), whereas class 4 to 6 should have low values - % (<0.15). -% old version ... -% weight(i) = min(1,max(0,2 * ( mean(res(i).wp(1:2)) - std(res(i).wp(1:2)) + ... -% mean(res(i).wp(4:6)) - std(res(i).wp(4:6))) + ... -% max(0,0.5 - 2 * sum( ... -% shiftdim(res(i).vr( res(i).lkp(:) < 7 )) ./ ... -% res(i).mn( res(i).lkp(:) < 7 )' .* ... -% res(i).mg( res(i).lkp(:) < 7 ) ) ) )); - weight(i) = min(1,max(0,1 - sum( shiftdim(res(i).vr) ./ ... - res(i).mn' .* res(i).mg ./ mean(res(i).mn(res(i).lkp(2)))) )); - - if numel(job.opts.tpm)>1 - fprintf('%s (fits%4.0f%%)',sprintf(repmat('\b',1,12)),weight(i)*100); - end - end - - - %% Mix of different TPMs - if numel(job.opts.tpm)>1 - % select most fitting TPM based on the affine registration - [tmp,weighti] = sort(weight,'descend'); clear tmp tpm; %#ok - - tpm{1} = spm_load_priors8(job.opts.tpm{weighti(1)}); - tpm{2} = spm_load_priors8(job.opts.tpm{weighti(2)}); - - if skullstripped - for i=1:2 - % update number of SPM gaussian classes - Ybg = 1 - spm_read_vols(tpm{i}.V(1)) - spm_read_vols(tpm{i}.V(2)) - spm_read_vols(tpm{i}.V(3)); - tpm{i}.V(4).dat = Ybg; - tpm{i}.dat{4} = Ybg; - tpm{i}.V(4).pinfo = repmat([1;0],1,size(Ybg,3)); - tpm{i}.dat(5:6) = []; - tpm{i}.V(5:6) = []; - tpm{i}.bg1(4) = tpm{i}.bg1(6); - tpm{i}.bg2(4) = tpm{i}.bg1(6); - tpm{i}.bg1(5:6) = []; - tpm{i}.bg2(5:6) = []; - end - end - - weightn = round(weight(weighti(1:2))/sum(weight(weighti(1:2)))*100)/100; - - if weightn(1)<0.9 - Affine = mAffine{end}; - stime = cat_io_cmd(sprintf(' Merge TPM %d and %d (%0.0f:%0.0f)',weighti(1:2),weightn*100),'g5','',job.extopts.verb-1,stime); - - obj.tpm = spm_load_priors8(job.opts.tpm{weighti(1)}); - for k1=1:min([numel(tpm{1}.dat),numel(tpm{2}.dat)]); - obj.tpm.dat{k1} = double(tpm{1}.dat{k1} * weightn(1) + ... - tpm{2}.dat{k1} * weightn(2)); - end - clear tpm1 tpm2; - - % final affine regisration with individual TPM - i = numel(job.opts.tpm) + 1; - mAffine{i} = spm_maff8(obj.image(1),... - obj.samp,obj.fwhm,obj.tpm,mAffine{weighti(1)},job.opts.affreg,40); - failedAffine = mean(abs(mAffine{end}(:) - mAffine{weighti(1)}(:))) < 0.2; - - % Final evaluation to test the result - if ~failedAffine - res(i).Affine = mAffine{end}; - Affine = mAffine{end}; - end - - res(i) = cat_spm_preproc8(obj); - weight(i) = min(1,max(0,1 - sum( shiftdim(res(i).vr) ./ ... - res(i).mn' .* res(i).mg ./ mean(res(i).mn(res(i).lkp(2)))) )); - fprintf('%s (fits%4.0f%%)',sprintf(repmat('\b',1,12)),weight(i)*100); - - if max(weight(1:end-1)) > weight(i) %failedAffine || - cat_io_cmd(sprintf(' Merging failed: Use only TPM %d!',weighti(1)),... - 'g5','',job.extopts.verb-1,stime); - - Affine = mAffine{weighti(1)}; - obj.tpm = spm_load_priors8(job.opts.tpm{weighti(1)}); - else - cat_io_cmd(' ','g5','',job.extopts.verb-1,stime); - end - else - obj.tpm = job.opts.tpm{weighti(1)}; - Affine = mAffine{weighti(1)}; - cat_io_cmd(sprintf(' Use only TPM %d.',weighti(1)),'g5','',job.extopts.verb-1,stime); - end - - - cat_io_cprintf('g5',sprintf('%5.0fs\n',etime(clock,stime0))); - cat_io_cmd(' ','','',job.extopts.verb-1); - - else - Affine = mAffine{1}; - end - - - tpm = obj.tpm; - -end \ No newline at end of file diff --git a/cat_tst_cattest.m b/cat_tst_cattest.m index c95e5aea..51599004 100644 --- a/cat_tst_cattest.m +++ b/cat_tst_cattest.m @@ -158,9 +158,7 @@ 'cat12_101_MAIN_segment' 1 1 'spm.tools.cat.estwrite.opts.tpm' 'tpm' ... {{fullfile(spm('dir'),'tpm','TPM.nii')} ... {fullfile(cat_get_defaults('extopts.pth_templates'),'TPM_Age11.5.nii')} ... % children - {fullfile(spm('dir'),'tpm','TPM.nii'); ... % multi-TPM - fullfile(cat_get_defaults('extopts.pth_templates'),'TPM_ba0_ha0.nii'); ... - fullfile(cat_get_defaults('extopts.pth_templates'),'TPM_ba0_ha0.nii')} }; + }; ... SPM bias ... SPM acc ... == parameter in default GUI == @@ -203,9 +201,7 @@ 'cat12_105_MAIN_segment_1173plus' 1 1 'spm.tools.cat.estwrite1173plus.opts.tpm' 'tpm' ... {{fullfile(spm('dir'),'tpm','TPM.nii')} ... {fullfile(cat_get_defaults('extopts.pth_templates'),'TPM_Age11.5.nii')} ... % children - {fullfile(spm('dir'),'tpm','TPM.nii'); ... % multi-TPM - fullfile(cat_get_defaults('extopts.pth_templates'),'TPM_ba0_ha0.nii'); ... - fullfile(cat_get_defaults('extopts.pth_templates'),'TPM_ba0_ha0.nii')} }; + }; ... SPM bias ... SPM acc ... == parameter in default GUI == diff --git a/cat_vol_pbt2.m b/cat_vol_pbt2.m deleted file mode 100644 index 408d79db..00000000 --- a/cat_vol_pbt2.m +++ /dev/null @@ -1,491 +0,0 @@ -function [Ygmt,Ypp,Ymf,Ywmd,Ywmdc] = cat_vol_pbt2(Ymf,opt) -% ______________________________________________________________________ -% -% Cortical thickness and surface position estimation. -% -% [Ygmt,Ypp,Ywmd,Ycsfd] = cat_vol_pbt(Ymf,opt) -% -% Ygmt: GM thickness map -% Ypp: percentage position map -% Ywmd: WM distance map -% Ycsfd: CSF distance map -% -% Ymf: tissue segment image or better the noise, bias, and -% intensity corrected -% -% opt.resV voxel resolution (only isotropic) -% opt.method choose of method {'pbt2x','pbt2'} with default=pbt2x as -% the method that is described in the paper. -% opt.pbtlas GM intensity correction to reduce myelination effects -% (added 201908) -% ______________________________________________________________________ -% -% Dahnke, R; Yotter R; Gaser C. -% Cortical thickness and central surface estimation. -% NeuroImage 65 (2013) 226-248. -% ______________________________________________________________________ -% -% Christian Gaser, Robert Dahnke -% Structural Brain Mapping Group (https://neuro-jena.github.io) -% Departments of Neurology and Psychiatry -% Jena University Hospital -% ______________________________________________________________________ -% $Id$ - - -% default variables and check/set function - if ~exist('opt','var'), opt=struct(); end - - def.resV = 1; - def.dmethod = 'eidist'; - def.method = 'pbt2x'; % pbt is worse ... just for tests! - def.debug = cat_get_defaults('extopts.verb')>2; - def.verb = cat_get_defaults('extopts.verb')-1; - def.pbtlas = 1; % no/yes - opt = cat_io_checkinopt(opt,def); - opt.resV = mean(opt.resV); - - minfdist = 2; - - if 0 - % RD 201803: - % Remove blood vessels & meninges - % This block is too aggressive and removes gyral peaks (thickness overestimation). - % However, some correction is required ... - Ymx = cat_vol_morph(Ymf>2.5,'l'); Ymf(Ymf>2.5 & ~Ymx)=2; clear Ymx - Ymx = cat_vol_morph(Ymf>2.2,'lo'); Ymf(Ymf>2.2 & ~Ymx)=2.1; clear Ymx - Ymx = cat_vol_morph(Ymf>1.5,'lo'); Ymf(Ymf>1.5 & ~Ymx)=1.2; clear Ymx - end - - - %% Distance maps - if (sum(round(Ymf(:))==Ymf(:)) / numel(Ymf))>0.9, bin=1; else bin=0; end - - % WM distance - % Estimate WM distance Ywmd and the outer CSF distance Ycsfdc to correct - % the values in CSF area are to limit the Ywmd to the maximum value that - % is possible within the cortex. - % The increasement of this area allow a more accurate and robust projection. - % cat_vol_eidist used speed map to align voxel to the closer gyrus - % that is not required for the correction map. - % - % RD 201803: - % The speed map weighting "max(0.5,min(1,Ymf/2))" is not strong enough to - % support asymmetric structures. The map "max(eps,min(1,((Ymf-1)/1.1).^4))" - % works much better but it leads to much higher thickness results (eg. in - % the Insula). - - newspeedmapF = 1; - - - if opt.verb, fprintf('\n'); end; stime2=clock; - YMM = cat_vol_morph(Ymf<1.5,'e',1) | isnan(Ymf); - switch opt.dmethod - case 'eidist' - % [D,I] = vbm_vol_eidist(B,L,[vx_vol,euclid,csf,setnan,verb]) - - if newspeedmapF - F = max(eps,min(1,((Ymf-1)/1.1).^4)); - else - F = max(0.5,min(1,Ymf/2)); % R1218 - end - - - if opt.pbtlas == 0 - % Simple pure intensity based model that show strong variation in - % myelinated regions. - YM = max(0,min(1,(Ymf-2))); YM(YMM) = nan; - stime = cat_io_cmd(' WM distance: ','g5','',opt.verb); - - else - % New local intensity optimization to remove myelination effects - % that works in a similar fasion as pbtlas but only focuses on the - % GM/WM boundary. - % This model was implemented after the results of Shahrzad Kharabian - % Masouleh (Juelich) that showed strong underestimations and increased - % variance in CAT compared to FS and CIVET. - % RD 201908 - - stime = cat_io_cmd(' WM boundar optimization: ','g5','',opt.verb); - vx_vol_org = 1; %max(1,sqrt(sum(opt.vmat(1:3,1:3).^2)))); - - - %% WM intensity: - % The result should be a map with an average value of 1 and no - % anatomical structures (maybe some differencence for regional - % folding independent pattern). - % Important to avoid gyral overestimations, however too small - % filter size / low number of iterations scan also lead to errors. - [Ymfr,resT2] = cat_vol_resize(Ymf,'reduceV',opt.resV,mean(vx_vol_org),32,'meanm'); - YM2max = max(0,min(1,Ymfr-2)); YM2max(YM2max<0.0) = 0; % lower threshold > higher filter range/more iterations - YM2max = YM2max .* cat_vol_morph(YM2max>0 ,'l',[10 0.1 ]); - for ix=1:3 %round(2/mean(resT2.vx_volr)) % iterations are faster/smoother than single calls with higher values - YM2max = cat_vol_localstat(YM2max,YM2max>0.0,min(10,2.0/mean(resT2.vx_volr)),3); % higher treshold and maximum filter - YM2max = cat_vol_localstat(YM2max,YM2max>0.0,min(10,1.0/mean(resT2.vx_volr)),1); % mean filter for smooth data - end - YM2max = cat_vol_approx(YM2max,'nh',resT2.vx_volr,2); % use (slow) high resolutions - YM2max = cat_vol_resize(YM2max,'dereduceV',resT2); - - % stime = cat_io_cmd(sprintf(' WM distance (GM - %0.2f): ',mean(vx_vol_org)),'g5','',opt.verb,stime); - %% GM intensity: - % find the higher intensity GM average and avoid the CSF PVE - lth = 1.8; %1.8; % 1.8 - 2.2 - [Ymfr,YM2maxr,resT2] = cat_vol_resize({Ymf,YM2max},'reduceV',opt.resV,mean(vx_vol_org),32,'meanm'); - - % This is a simple tissue thickness model from a more central GM layer - % (currently 2.1, possible range 1.5 to 2.2) to avoid corrections - %% in thicker areas (mostly gyris). - Ygd = cat_vbdist(2.1 - Ymfr); % estimate distance map to central/WM surface, lower thresholds are also possible (range 1.5 to 2.2, default 2.1) - Ygdt = cat_vol_pbtp2(max(2,min(3,4-Ymfr)),Ygd,inf(size(Ygd),'single')) * mean(resT2.vx_volr); - Ygdt = cat_vol_median3(Ygdt,Ygdt>0.01,Ygdt>0.01); - Ygdt = cat_vol_localstat(Ygdt,Ygdt>0.1,1/mean(resT2.vx_volr),1); - - Ywd = cat_vbdist(2.5 - Ymfr); - Ywdt = cat_vol_pbtp2(max(2,min(3,4-Ymfr)),Ywd,inf(size(Ywd),'single')) * mean(resT2.vx_volr); - Ywdt = cat_vol_median3(Ywdt,Ywdt>0.01,Ywdt>0.01); - Ywdt = cat_vol_localstat(Ywdt,Ywdt>0.1,1/mean(resT2.vx_volr),1); - - - %% - % estimate local GM average intensity - YM2min = max(0,min(1.5,Ymfr - lth)) .* ... - (Ymfr<( (YM2maxr+2) .* (2.5 + min(0.45,Ygdt/8) )) /3) .* ... - (Ymfr<( (YM2maxr+2) .* (2.5 + min(0.45,Ywdt/8) )) /3) .* ... - (Ywd<4*mean(resT2.vx_volr)) .*... - (cat_vbdist(single(Ymfr./YM2maxr<2))<3.5/mean(resT2.vx_volr)); %clear YM2maxr Ywdt - - % remove small dots - YM2min = YM2min .* cat_vol_morph(YM2min>0 ,'l',[100 0.01]); - YM2min(smooth3(YM2min>0)<0.4) = 0; - YM2min = cat_vol_median3(YM2min,YM2min>0,YM2min>0); - % avoid PVE - YM2min2 = cat_vol_localstat(YM2min,YM2min>0,1,2); - YM = Ymfr>2 & YM2min==0; YM = cat_vol_morph(YM,'dd',1) & ~YM; YM2min( YM ) = YM2min2( YM ); - YM2min2 = cat_vol_localstat(YM2min,YM2min>0,1,3); - YM = Ymfr<2 & YM2min==0; YM = cat_vol_morph(YM,'dd',1) & ~YM; YM2min( YM ) = YM2min2( YM ); - %clear Ymfr YM - - % remove further atypical values - YM2mina = YM2min; - for ix=1:10 - YM2mina = cat_vol_localstat(YM2mina,YM2mina>0,2/mean(resT2.vx_volr),1); - end - YM2min = YM2min .* ((YM2mina - YM2min)>-0.2 & (YM2mina - YM2min)<0.3); clear YM2mina; - - % remove PVE - YM2min = cat_vol_median3(YM2min,YM2min>0); - - % further local filtering - for ix=1:10 %round(5/mean(resT2.vx_volr)) - YM2min = cat_vol_localstat(YM2min,YM2min>0,1/mean(resT2.vx_volr),1); - end - YM2min(Ywd>6*mean(resT2.vx_volr)) = 0.01; - - YM2min = cat_vol_approx(YM2min,'nh',resT2.vx_volr,2); % use (slow) high resolutions - YM2min = cat_vol_resize(YM2min,'dereduceV',resT2); - Ywd = cat_vol_resize(Ywd,'dereduceV',resT2); - YM2min = YM2min - (2 - lth); - %stime = cat_io_cmd(sprintf(' WM distance (GM - %0.2f): ',mean(vx_vol_org)),'g5','',opt.verb,stime); - %% create intensity optimized boundary image - YM = max(0,min(1,((Ymf - (YM2min + 2)) ./ (YM2max - YM2min) ))); - %% use only fat areas - YMo = max(0,min(1,(Ymf-2))); - YMc = YMo - YM; - YMc(YMc>0 & YMc<0.15) = YMc(YMc>0 & YMc<0.15); - YMc(YMc>0 & YMc<0.10) = YMc(YMc>0 & YMc<0.10)/2; - YMc(YMc>0 & YMc<0.05) = YMc(YMc>0 & YMc<0.05)/4; - YMc(YMc<0.05) = 0; - YM = YMo - YMc; - clear YMo YMc YMcs; - - % distance corrections to stabilize thin structures - YM = max(YM,min(1,Ywd/mean(resT2.vx_volr)/2 - 2)); % L4 distance based corrections - YMd = 0.75 - YM; YMd(YM>0.9) = nan; - YM = max(YM,min(1,cat_vol_eidist(YMd,ones(size(YM),'single'),[1 1 1],1,1,0,opt.debug) - 0.5)); - %stime = cat_io_cmd(sprintf(' WM distance (GM - %0.2f): ',mean(vx_vol_org)),'g5','',opt.verb,stime); - %% - for ix=1:3, YM(YM>0.5 & smooth3(YM>0.5)<0.2) = 0.45; end - YM(cat_vol_morph(YM>0.9 ,'lc')) = 1; - YM(YMM) = nan; - clear YM2max YM2min - - stime = cat_io_cmd(' WM distance: ','g5','',opt.verb,stime); - end - - - %% final distance estimation - YMwm = YM; - Ywmd = cat_vol_eidist(YM,F,[1 1 1],1,1,0,opt.debug); - - if newspeedmapF - F = max(eps,min(1,((Ymf-1)/1.1).^4)); - else - F = max(1.0,min(1,Ymf/2)); % R1218 - no speed difference! - end - YM = max(0,min(1,(Ymf-1))); YM(YMM) = nan; Ywmdc = cat_vol_eidist(YM,F,[1 1 1],1,1,0,opt.debug); - - clear F; - - case 'vbdist' - stime = cat_io_cmd(' WM distance: ','g5','',opt.verb); stime2=stime; - YM = max(0,min(1,(Ymf-2))); Ywmd = max(0,cat_vbdist(single(YM>0.5),~YMM)-0.5); - YM = max(0,min(1,(Ymf-1))); Ywmdc = max(0,cat_vbdist(single(YM>0.5),~YMM)-0.5); - end - clear YMM; - if ~bin - % limit the distance values outside the GM/CSF boudary to the distance possible in the GM - YM = Ywmd>minfdist & Ymf<=1.5; Ywmd(YM) = Ywmd(YM) - Ywmdc(YM); Ywmd(isinf(Ywmd)) = 0; clear Ycsfdc; - % smoothing of distance values inside the GM - %YM = Ywmd>minfdist & Ymf> 1.5; YwmdM = Ywmd; YwmdM = cat_vol_localstat(YwmdM,YM,1,1); Ywmd(YM) = YwmdM(YM); - % smoothing of distance values outside the GM - YM = Ywmd>minfdist & Ymf<=1.5; YwmdM = Ywmd; for i=1:1, YwmdM = cat_vol_localstat(YwmdM,YM,1,1); end; Ywmd(YM) = YwmdM(YM); - % reducing outliers in the GM/CSF area - YM = Ywmd>minfdist & Ymf< 2.0; YwmdM = Ywmd; YwmdM = cat_vol_median3(YwmdM,YM,YM); Ywmd(YM) = YwmdM(YM); clear YwmdM YM; - end - - minfdist = 2/opt.resV; - % CSF distance - % Similar to the WM distance, but keep in mind that this map is - % incorrect in blurred sulci that is handled by PBT - stime = cat_io_cmd(' CSF distance: ','g5','',opt.verb,stime); - if exist('YMwm','var') - YMM = cat_vol_morph(Ymf<1.5,'e',1) | cat_vol_morph(YMwm>0.5,'e',1) | isnan(Ymf); % this was dilate??? - else - YMM = cat_vol_morph(Ymf<1.5,'e',1) | cat_vol_morph(Ymf>2.5,'e',1) | isnan(Ymf); % this was dilate??? - end - switch opt.dmethod - case 'eidist' - if newspeedmapF - F = max(eps,min(1,(4-Ymf)/2).^2); - else - F = max(0.5,min(1,(4-Ymf)/2)); % R1218 - end - YM = max(0,min(1,(2-Ymf))); YM(YMM) = nan; Ycsfd = cat_vol_eidist(YM,F,[1 1 1],1,1,0,opt.debug); - if newspeedmapF - F = max(eps,min(1,(4-Ymf)/2).^2); - else - F = max(1,min(1,(4-Ymf)/2)); % R1218 - end - if exist('YMwm','var') - YM = 1 - YMwm; YM(YMM) = nan; Ywmdc = cat_vol_eidist(YM,F,[1 1 1],1,1,0,opt.debug); Ywmdx = Ywmdc; - else - YM = max(0,min(1,(3-Ymf))); YM(YMM) = nan; Ywmdc = cat_vol_eidist(YM,F,[1 1 1],1,1,0,opt.debug); - YM = max(0,min(1,(2.7-Ymf))); YM(YMM) = nan; Ywmdx = cat_vol_eidist(YM,F,[1 1 1],1,1,0,opt.debug)+0.3; - end - clear F; - case 'vbdist' - YM = max(0,min(1,(2-Ymf))); Ycsfd = max(0,cat_vbdist(single(YM>0.5),~YMM)-0.5); - YM = max(0,min(1,(3-Ymf))); Ywmdc = max(0,cat_vbdist(single(YM>0.5),~YMM)-0.5); - YM = max(0,min(1,(2.7-Ymf))); Ywmdx = max(0,cat_vbdist(single(YM>0.5),~YMM)-0.2); - end - Ywmdc = min(Ywmdc,Ywmdx); - clear YMM Ywmdx; - if ~bin - if exist('YMwm','var') - YM = Ycsfd>minfdist & YMwm>=0.5; Ycsfd(YM) = Ycsfd(YM) - Ywmdc(YM); Ycsfd(isinf(-Ycsfd)) = 0; clear Ywmdc; - %YM = Ycsfd>minfdist & YMwm< 0.5; YcsfdM = Ycsfd; YcsfdM = cat_vol_localstat(YcsfdM,YM,1,1); Ycsfd(YM) = YcsfdM(YM); - YM = Ycsfd>minfdist & YMwm>=0.5; YcsfdM = Ycsfd; for i=1:1, YcsfdM = cat_vol_localstat(YcsfdM,YM,1,1); end; Ycsfd(YM) = YcsfdM(YM); - else - YM = Ycsfd>minfdist & Ymf>=2.5; Ycsfd(YM) = Ycsfd(YM) - Ywmdc(YM); Ycsfd(isinf(-Ycsfd)) = 0; clear Ywmdc; - %YM = Ycsfd>minfdist & Ymf< 2.5; YcsfdM = Ycsfd; YcsfdM = cat_vol_localstat(YcsfdM,YM,1,1); Ycsfd(YM) = YcsfdM(YM); - YM = Ycsfd>minfdist & Ymf>=2.5; YcsfdM = Ycsfd; for i=1:1, YcsfdM = cat_vol_localstat(YcsfdM,YM,1,1); end; Ycsfd(YM) = YcsfdM(YM); - end - YM = Ycsfd>minfdist & Ymf> 2.0; YcsfdM = Ycsfd; YcsfdM = cat_vol_median3(YcsfdM,YM,YM); Ycsfd(YM) = YcsfdM(YM); clear YcsfdM YM; - end - - - %% PBT thickness mapping - % PBT is the default thickness estimation, but PBT2x is the optimized - % version that usex both sulci and gyri refinements, because not only - % thin sulci can be blurred. PBT2x is also the method that is - % described in the paper. - % PBTv is new version that uses the volume rather than the distance. - % Although this works in principle, this is biased by interpolation - % artifacts and not-optimal WMD mapping. - iter = 0;%round(0.5/mean(opt.resV)); - fs = 0.5; - if strcmp(opt.method,'pbtv') - Ywmdo = Ywmd+0; - - %% estimate cortical thickness and map the local volumes - stime = cat_io_cmd(' PBTV thickness: ','g5','',opt.verb,stime); - [Ygmt,Yv1,Yv2,Ypp] = cat_vol_pbtv(Ymf,Ywmd,Ycsfd); - Ygmts = Ygmt; for i=1:iter, Ygmts = cat_vol_localstat(Ygmts,(Ygmt>1 | Ypp>0.1) & Ygmt>0 & (Ygmt>1 | Ymf>1.8),1,1); end; Ygmt(Ygmts>0) = Ygmts(Ygmts>0); - - % volume-based PP map - Ypp = ((Ygmt - Ywmd)./(Ygmt + eps) + (Ymf>2.5))*0.5 + 0.5*min( 1, Yv1./(Yv1 + Yv2 + eps) + (Ymf>2.5)); - - - elseif strcmp(opt.method,'pbt2x') || strcmp(opt.method,'pbt2x2') - % Estimation of the cortical thickness with sulcus (Ygmt1) and gyri - % correction (Ygmt2) to create the final thickness as the minimum map - % of both. - stime = cat_io_cmd(' PBT2x thickness: ','g5','',opt.verb,stime); - - % estimate thickness with PBT approach - if opt.pbtlas, Ymfo=Ymf; Ymf = single(1 + 2*((Ywmd>0 & Ycsfd>0) | Ymfo>2)) - (Ywmd>0 & Ycsfd>0); end - Ygmt1 = cat_vol_pbtp2((Ymf),Ywmd,Ycsfd); - Ygmt2 = cat_vol_pbtp2((4-Ymf),Ycsfd,Ywmd); - - % Error handling - % For some unkown reasons the sulcus reconstruction of cat_vol_pbtp failed in some cases (not directly reproducable). - % Reprocessing is solving this problem, but further investigation of cat_vol_pbtp.cpp would be good (RD 20190811). - % Maybe it depends on the initialization of the regions, e.g., using Ymf without rounding and incorrect boundary seems to increase the problems. - % Now use cat_vol_pbtp2.cpp (RD20200111). - mask = @(Y) Y(:)>0 & Y(:)<1000000; - rerun = 0; rerunlim = 3; - while rerun <= rerunlim && isnan( mean( Ygmt1(mask(Ygmt1))) ) || mean( Ygmt1(mask(Ygmt1)))>100 - Ygmt1 = cat_vol_pbtp2((Ymf),Ywmd,Ycsfd); - rerun = rerun + 1; - pause(rand*3); - end - if rerun == rerunlim && isnan( mean( Ygmt1(mask(Ygmt1))) ) || mean( Ygmt1(mask(Ygmt1)))>100 - error('cat_vol_pbtp2:bad_mapping1','Untypcial values in PBT thickness mapping detected. '); - end - rerun = 0; - while rerun <= rerunlim && isnan( mean( Ygmt2(mask(Ygmt2))) ) || mean( Ygmt2(mask(Ygmt2)))>100 - Ygmt2 = cat_vol_pbtp2((4-Ymf),Ycsfd,Ywmd); - rerun = rerun + 1; - pause(rand*3); - end - if rerun == rerunlim && rerunlim && isnan( mean( Ygmt2(mask(Ygmt2))) ) || mean( Ygmt2(mask(Ygmt2)))>100 - error('cat_vol_pbtp2:bad_mapping2','Untypcial values in PBT thickness mapping detected. '); - end - - %% avoid meninges ! - update_WMD = strcmp(opt.method,'pbt2x2'); - - Ygmt1 = min(Ygmt1,Ycsfd+Ywmd); - if update_WMD % eigentlich muss das hier an sein ... - Ygmt2 = min(Ygmt2 + (Ygmt2>0) .* 2,Ycsfd+Ywmd); - else - Ygmt2 = min(Ygmt2,Ycsfd+Ywmd); - end - - % median filter to remove extreme outliers - YM = Ymf>1.5 & Ymf<2.5; - Ygmt1 = cat_vol_median3(Ygmt1,YM,Ygmt1>0,0.2); - Ygmt2 = cat_vol_median3(Ygmt2,YM,Ygmt2>0,0.2); - - %% estimation of Ypp for further GM filtering without sulcal blurring - if update_WMD - [Ygmt,Yi] = min(cat(4,Ygmt1,Ygmt2),[],4); - Ywmdc = (Ywmd.*(Yi==1) + (Yi==2).*(Ygmt2 - Ycsfd)); - Ywmdc = cat_vol_median3(Ywmdc,Yi==2,Yi>0,0.05); - Ywmdcs = cat_vol_localstat(Ywmdc,Yi==2,2,1); Ywmdc(Ywmdcs>0) = Ywmdc(Ywmdcs>0); clear Ycsfdcs - Ywmd = Ywmdc; - else - Ygmt = min(cat(4,Ygmt1,Ygmt2),[],4); - end - - Ygmt = Ygmt*(1-fs) + fs*cat_vol_median3(Ygmt,Ygmt>0,Ygmt>0,0.05); - Ypp = zeros(size(Ymf),'single'); Ypp(Ymf>=2.5) = 1; YM = Ygmt>0; - Ypp(YM) = min(Ycsfd(YM),Ygmt(YM) - Ywmd(YM)) ./ (Ygmt(YM) + eps); Ypp(Ypp>2) = 0; - Ygmts = cat_vol_localstat(Ygmt,Ygmt>0 & Ypp>0.1,1,1); Ygmt(Ygmts>0) = Ygmt(Ygmts>0)*(1-fs) + fs*Ygmts(Ygmts>0); clear Ygmtps - %% - Ypp(Ypp==0 & Ymf>=2.5 & Ywmd<=0) = 1; - Ypp(Ypp==0 & Ymf>=2.5 & Ywmd<1) = 0.99; - Ypp(Ypp==0 & Ymf>=2.0 & Ywmd<2) = 0.985; - Ypp(Ypp==0 & Ymf>=2.5) = 0.95; - - Ypp = Ypp*(1-fs) + fs*cat_vol_median3(Ypp,Ymf>1.1 & Ymf<2.9 & (Ypp>0.1 | Ypp==0),true(size(Ypp)),0.1); - Ypp = Ypp*(1-fs) + fs*cat_vol_median3(Ypp,Ymf>1.1 & Ymf<2.9 & (Ypp>0.1 | Ypp==0),true(size(Ypp)),0.1); - - Ypp(Ypp>0.5 & ~cat_vol_morph(Ypp>0.5,'ldo',0.50/opt.resV)) = 0.49; - Ypp(Ypp<0.5 & ~cat_vol_morph(Ypp<0.5,'l' ,0.25/opt.resV)) = 0.51; - - - else - % Estimation of thickness map Ygmt and percentual position map Ypp. - stime = cat_io_cmd(' PBT2 thickness: ','g5','',opt.verb,stime); - - % estimate thickness with PBT approach - if opt.pbtlas, Ymfo=Ymf; Ymf = single(1 + 2*((Ywmd>0 & Ycsfd>0) | Ymfo>2)) - (Ywmd>0 & Ycsfd>0); end - - % estimate thickness with PBT approach - if opt.pbtlas, Ymfo=Ymf; Ymf = single(1 + 2*((Ywmd>0 & Ycsfd>0) | Ymfo>2)) - (Ywmd>0 & Ycsfd>0); end - Ygmt = cat_vol_pbtp2( (Ymf) ,Ywmd,Ycsfd); - Ygmt = Ygmt*(1-fs) + fs*cat_vol_median3(Ygmt,Ygmt>0,Ygmt>0,0.05); - Ypp = zeros(size(Ymf),'single'); Ypp(Ymf>=2.5) = 1; YM = Ygmt>0; - Ypp(YM) = min(Ycsfd(YM),Ygmt(YM) - Ywmd(YM)) ./ (Ygmt(YM) + eps); Ypp(Ypp>2) = 0; - Ygmts = cat_vol_localstat(Ygmt,Ygmt>0 & Ypp>0.1,1,1); Ygmt(Ygmts>0) = Ygmt(Ygmts>0)*(1-fs) + fs*Ygmts(Ygmts>0); clear Ygmtps - - % Error handling - % For some unkown reasons the sulcus reconstruction of cat_vol_pbtp failed in some cases (not directly reproducable). - % Reprocessing is solving this problem, but further investigation of cat_vol_pbtp.cpp would be good (RD 20190811). - % Maybe it depends on the initialization of the regions, e.g., using Ymf without rounding and incorrect boundary seems to increase the problems. - % Now use cat_vol_pbtp2.cpp (RD20200111). - mask = @(Y) Y(:)>0 & Y(:)<1000000; - rerun = 0; rerunlim = 3; - while rerun <= rerunlim && isnan( mean( Ygmt(mask(Ygmt))) ) || mean( Ygmt(mask(Ygmt)))>100 - Ygmt = cat_vol_pbtp2((Ymf),Ywmd,Ycsfd); - rerun = rerun + 1; - pause(rand*3); - end - if rerun == rerunlim && isnan( mean( Ygmt(mask(Ygmt))) ) || mean( Ygmt(mask(Ygmt)))>100 - error('cat_vol_pbtp2:bad_mapping','Untypcial values in PBT thickness mapping detected. '); - end - - % avoid meninges ! - Ygmt = min(Ygmt,Ycsfd+Ywmd); - - % median filter to remove outliers - Ygmt = cat_vol_median3(Ygmt,Ygmt>0,Ygmt>0); - - %YM = Ymf>=1.5 & Ymf<2.5; Ypp = zeros(size(Ymf),'single'); Ypp(Ymf>=2.5) = 1; - %Ypp(YM) = min(Ycsfd(YM),Ygmt(YM) - Ywmd(YM)) ./ (Ygmt(YM) + eps); Ypp(Ypp>2) = 0; - - Ypp(Ypp==0 & Ymf>=2.5 & Ywmd<=0) = 1; - Ypp(Ypp==0 & Ymf>=2.5 & Ywmd<1) = 0.99; - Ypp(Ypp==0 & Ymf>=2.0 & Ywmd<2) = 0.985; - Ypp(Ypp==0 & Ymf>=2.5) = 0.95; - - Ypp = Ypp*(1-fs) + fs*cat_vol_median3(Ypp,Ymf>1.1 & Ymf<2.9 & (Ypp>0.1 | Ypp==0),true(size(Ypp)),0.1); - Ypp = Ypp*(1-fs) + fs*cat_vol_median3(Ypp,Ymf>1.1 & Ymf<2.9 & (Ypp>0.1 | Ypp==0),true(size(Ypp)),0.1); - - Ypp(Ypp>0.5 & ~cat_vol_morph(Ypp>0.5,'ldo',0.50/opt.resV)) = 0.49; - Ypp(Ypp<0.5 & ~cat_vol_morph(Ypp<0.5,'l' ,0.25/opt.resV)) = 0.51; - - - % filter result - if exist('Ymfo','var'); Ymf=Ymfo; end - Ygmts = Ygmt; for i=1:iter, Ygmts = cat_vol_localstat(Ygmts,(Ygmt>1 | Ypp>0.1) & Ygmt>0 & (Ygmt>1 | Ymf>1.8),1,1); end; Ygmt(Ygmts>0) = Ygmts(Ygmts>0); - - % filter result - Ygmts = Ygmt; for i=1:iter, Ygmts = cat_vol_localstat(Ygmts,(Ygmt>1 | Ypp>0.1) & Ygmts>0 & (Ygmt>1 | Ymf>1.8),1,1); end; Ygmt(Ygmts>0) = Ygmts(Ygmts>0); - end - - - - - - - %% Estimation of a mixed percentual possion map Ypp. - if nargout>3 || ~strcmp(opt.method,'pbtv') - YM = Ymf>=1.5 & Ymf<2.5 & Ygmt>eps; - Ywmdc = Ycsfd; Ywmdc(YM) = min(Ycsfd(YM),Ygmt(YM) - Ywmd(YM)); - end - if ~strcmp(opt.method,'pbtv') - Ypp = zeros(size(Ymf),'single'); Ypp(Ymf>=2.5)=1; - Ypp(YM) = Ywmdc(YM) ./ (Ygmt(YM) + eps); - Ypp(Ypp>2) = 0; - end - - % YM = (Ygmt<=opt.resV & Ywmd<=opt.resV & Ygmt>0); Ypp(YM) = (Ymf(YM)-1)/2 - 0.2; % correction of voxel with thickness below voxel resolution - Ypp(isnan(Ypp)) = 0; - Ypp(Ypp<0) = 0; - - %% Final corrections for position map with removing of non brain objects. - % ds('d2','',1,Ymf/3,Ywmd/3,Ygmt/5,Ypp,70) - - % Final corrections for thickness map with thickness limit of 10 mm. - % Resolution correction of the thickness map after all other operations, - % because PBT actually works only with the voxel-distance (isotropic 1 mm) - Ygmt = Ygmt*opt.resV; - Ygmt(Ygmt>10) = 10; - if exist('Ycsfdc','var'), Ywmdc = Ywmdc*opt.resV; end - if exist('Ywmd','var'), Ywmd = Ywmd*opt.resV; end - - cat_io_cmd(' ','g5','',opt.verb,stime); - if opt.debug, cat_io_cmd(' ','','',opt.debug,stime2); end -end \ No newline at end of file diff --git a/cat_vol_pbtp.c b/cat_vol_pbtp.c index 9706c6e8..cb5a362c 100644 --- a/cat_vol_pbtp.c +++ b/cat_vol_pbtp.c @@ -209,7 +209,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { } } - for (int i=0;iopt.HB) GMT[i]=0; + for (int i=0;iopt.HB) GMT[i] = 0.0; @@ -217,7 +217,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* final settings... ======================================================================= */ - float CSFDc = 0, GMTi, CSFDi; /* 0.125 */ + float CSFDc = 0.0, GMTi, CSFDi; /* 0.125 */ for (int i=0;i=opt.LB && SEG[i]<=opt.LB) { diff --git a/cat_vol_pbtp2.c b/cat_vol_pbtp2.c deleted file mode 100644 index 5ceed40a..00000000 --- a/cat_vol_pbtp2.c +++ /dev/null @@ -1,248 +0,0 @@ -/* quasi-euclidean distance calculation - * _____________________________________________________________________________ - * [GMT,RPM,WMD,CSFD,II] = cat_vol_pbtp(SEG,WMD,CSFD[,opt]) - * - * SEG = (single) segment image with low and high boundary bd - * GMT = (single) thickness image - * RPM = (single) radial position map - * WMD = (single) CSF distance map - * CSFD = (single) CSF distance map - * II = (uint32) index of the inner (WM) boundary voxel - * - * opt.bd = (single) [low,high] boundary values (default 1.5 and 2.5) - * opt.CSFD = calculate CSFD - * opt.PVE = use PVE information (0=none,1=fast,2=exact) - * - * TODO: - * - eikonal distance for subsegmentation (region growing) - * - own labeling ( - * ______________________________________________________________________ - * - * Christian Gaser, Robert Dahnke - * Structural Brain Mapping Group (https://neuro-jena.github.io) - * Departments of Neurology and Psychiatry - * Jena University Hospital - * ______________________________________________________________________ - * $Id$ - */ - -#include "mex.h" -#include "math.h" -#include - -struct opt_type { - int CSFD; /* use CSFD */ - int PVE; /* 0, 1=fast, 2=exact */ - float LB, HB, LLB, HLB, LHB, HHB; /* boundary */ - int sL[3]; - // ... - } opt; - -float min(float a, float b) { - if (ab) return a; else return b; -} - -// get all values of the voxels which are in WMD-range (children of this voxel) -float pmax(const float GMT[], const float RPM[], const float SEG[], const float ND[], const float WMD, const float SEGI, const int sA) { - float T[27]; - float n=0.0, maximum=WMD; - - for (int i=0;i<27;i++) T[i]=-1; - - /* the pure maximum */ - /* (GMT[i]<1e15) && (maximum < GMT[i]) && ((RPM[i]-ND[i]*1.25)<=WMD) && ((RPM[i]-ND[i]*0.5)>WMD) && (SEGI)>=SEG[i] && SEG[i]>1 && SEGI>1.66) */ - for (int i=0;i<=sA;i++) { - if ( ( GMT[i] < 1e15 ) && ( maximum < GMT[i] ) && /* thickness/WMD of neighbors should be larger */ - ( SEG[i] >= 1.0 ) && ( SEGI>1.5 && SEGI<=2.5 ) && /* projection range */ - ( ( ( RPM[i] - ND[i] * min( 1.4, 1.1 + max(0,SEGI - SEG[i])) ) <= WMD ) || ( SEG[i]<1.5 ) ) && /* upper boundary - maximum distance - be more tollerant with the distance if the intensiy is lower */ - ( ( ( RPM[i] - ND[i] * min( 0.3, 0.6 - max(0,SEGI - SEG[i])) ) > WMD ) || ( SEG[i]<1.5 ) ) && /* lower boundary - minimum distance - corrected values outside */ - ( ( ( (SEGI * max(1.01,min(1.1,SEGI-1.5)) ) >= SEG[i] ) ) || ( SEG[i]<1.5 ) ) ) /* for high values will project data over sulcal gaps */ - { maximum = GMT[i]; } - } - - - /* the mean of the highest values*/ - float maximum2=maximum; float m2n=0.0; - for (int i=0;i<=sA;i++) { - if ( ( GMT[i] < 1e15 ) && ( maximum < GMT[i] ) && - ( SEG[i] >= 1.0 ) && ( SEGI>1.5 && SEGI<=2.5 ) && - ( ( ( RPM[i] - ND[i] * min( 1.4, 1.1 + max(0,SEGI - SEG[i])) ) <= WMD ) || ( SEG[i]<1.5 ) ) && /* upper boundary - maximum distance - be more tollerant with the distance if the intensiy is lower */ - ( ( ( RPM[i] - ND[i] * min( 0.3, 0.6 - max(0,SEGI - SEG[i])) ) > WMD ) || ( SEG[i]<1.5 ) ) && /* lower boundary - minimum distance - corrected values outside */ - ( ( ( (SEGI * max(1.01,min(1.1,SEGI-1.5)) ) >= SEG[i] ) ) || ( SEG[i]<1.5 ) ) ) - { maximum2 = maximum2 + GMT[i]; m2n++; } - } - if ( m2n > 0.0 ) maximum = (maximum2 - maximum)/m2n; - - return maximum; -} - - - - -/* estimate x,y,z position of index i in an array size sx,sxy=sx*sy... */ -void ind2sub(int i, int *x, int *y, int *z, int snL, int sxy, int sy) { - /* not here ... - * if (i<0) i=0; - * if (i>=snL) i=snL-1; - */ - - *z = (int)floor( (double)i / (double)sxy ) ; - i = i % (sxy); - *y = (int)floor( (double)i / (double)sy ) ; - *x = i % sy ; -} - - -// main function -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs<3) mexErrMsgTxt("ERROR: not enough input elements\n"); - if (nrhs>4) mexErrMsgTxt("ERROR: too many input elements.\n"); - if (nlhs<2) mexErrMsgTxt("ERROR: not enough output elements.\n"); - if (nlhs>2) mexErrMsgTxt("ERROR: too many output elements.\n"); - if (mxIsSingle(prhs[0])==0) mexErrMsgTxt("ERROR: first input must be an 3d single matrix\n"); - - - // main information about input data (size, dimensions, ...) - const mwSize *sL = mxGetDimensions(prhs[0]); - mwSize sSEG[] = {sL[0],sL[1],sL[2]}; - const int dL = mxGetNumberOfDimensions(prhs[0]); - const int nL = mxGetNumberOfElements(prhs[0]); - const int x = sL[0]; - const int y = sL[1]; - const int xy = x*y; - const float s2 = sqrt(2.0); - const float s3 = sqrt(3.0); - const int nr = nrhs; - - // indices of the neighbor Ni (index distance) and euclidean distance NW - const int NI[] = { 0, -1,-x+1, -x,-x-1, -xy+1,-xy,-xy-1, -xy+x+1,-xy+x,-xy+x-1, -xy-x+1,-xy-x,-xy-x-1}; - const float ND[] = {0.0,1.0, s2,1.0, s2, s2,1.0, s2, s3, s2, s3, s3, s2, s3}; - const int sN = sizeof(NI)/4; - float DN[sN],DI[sN],GMTN[sN],WMDN[sN],SEGN[sN],DNm; - - float du, dv, dw, dnu, dnv, dnw, d, dcf, WMu, WMv, WMw, GMu, GMv, GMw, SEGl, SEGu, tmpfloat; - int ni,u,v,w,nu,nv,nw, tmpint, WMC=0, CSFC=0; - - // main volumes - actual without memory optimization ... - plhs[0] = mxCreateNumericArray(dL,sL,mxSINGLE_CLASS,mxREAL); - plhs[1] = mxCreateNumericArray(dL,sL,mxSINGLE_CLASS,mxREAL); - -/* not yet defined - plhs[2] = mxCreateNumericArray(dL,sL,mxSINGLE_CLASS,mxREAL); - plhs[3] = mxCreateNumericArray(dL,sL,mxSINGLE_CLASS,mxREAL); - plhs[4] = mxCreateNumericArray(dL,sL,mxUINT32_CLASS,mxREAL); -*/ - - // input variables - float*SEG = (float *)mxGetPr(prhs[0]); - float*WMD = (float *)mxGetPr(prhs[1]); - float*CSFD = (float *)mxGetPr(prhs[2]); - - /*if ( nrhs>1) { - tmpint = (int)mxGetScalar(mxGetField(prhs[1],1,"CSFD")); printf("X=%d", tmpint); if ( tmpint!=NULL && (tmpint>=0 && tmpint<=1) ) opt.CSFD = tmpint; else opt.CSFD = 1; - tmpint = (int)mxGetScalar(mxGetField(prhs[1],1,"PVE")); printf("X=%d", tmpint); if ( tmpint!=NULL && (tmpint>=0 && tmpint<=2) ) opt.PVE = tmpint; else opt.PVE = 2; - tmpfloat = (float)mxGetScalar(mxGetField(prhs[1],1,"LB")); printf("X=%d", tmpfloat); if ( tmpfloat!=NULL ) opt.LB = tmpfloat; else opt.LB = 1.5; - tmpfloat = (float)mxGetScalar(mxGetField(prhs[1],1,"HB")); printf("X=%d", tmpfloat); if ( tmpfloat!=NULL ) opt.HB = tmpfloat; else opt.HB = 2.5; - } - else */{ opt.CSFD = 1;opt.PVE = 2;opt.LB = 1.5;opt.HB = 2.5; } - opt.LLB=floor(opt.LB), opt.HLB=ceil(opt.LB), opt.LHB=floor(opt.HB), opt.HHB=ceil(opt.HB); - - // output variables - float *GMT = (float *)mxGetPr(plhs[0]); - float *RPM = (float *)mxGetPr(plhs[1]); - - // intitialisiation - for (int i=0;i=opt.HB ) WMC++; - if ( SEG[i]<=opt.LB ) CSFC++; - } - if (WMC==0) mexErrMsgTxt("ERROR: no WM voxel\n"); - if (CSFC==0) opt.CSFD = 0; - - - -// thickness calcuation -// ============================================================================= - for (int i=0;iopt.LLB && SEG[i]=nL) || (abs(nu-u)>1) || (abs(nv-v)>1) || (abs(nw-w)>1)) ni=i; - GMTN[n] = GMT[ni]; WMDN[n] = RPM[ni]; SEGN[n] = SEG[ni]; - } - - // find minimum distance within the neighborhood - DNm = pmax(GMTN,WMDN,SEGN,ND,WMD[i],SEG[i],sN); - GMT[i] = DNm; - } - } - - for (int i=nL-1;i>=0;i--) { - if (SEG[i]>opt.LLB && SEG[i]=nL) || (abs(nu-u)>1) || (abs(nv-v)>1) || (abs(nw-w)>1)) ni=i; - GMTN[n] = GMT[ni]; WMDN[n] = RPM[ni]; SEGN[n] = SEG[ni]; - } - - // find minimum distance within the neighborhood - DNm = pmax(GMTN,WMDN,SEGN,ND,WMD[i],SEG[i],sN); - if ( GMT[i] < DNm && DNm>0 ) GMT[i] = DNm; - } - } - - for (int i=0;iopt.HB) GMT[i]=0.0; //WMD[i] - - - - - -// final setings... -// ============================================================================= - float CSFDc = 0.0, GMTi, CSFDi; // 0.125 - for (int i=0;i=opt.LB & SEG[i]<=opt.LB) { - GMTi = CSFD[i] + WMD[i]; - CSFDi = GMT[i] - WMD[i]; - - if ( CSFD[i]>CSFDi ) CSFD[i] = CSFDi; - else GMT[i] = GMTi; - } - } - - -// estimate RPM -// ============================================================================= - for (int i=0;i=opt.HB ) - RPM[i]=1.0; - else { - if ( SEG[i]<=opt.LB || GMT[i]==0.0 ) - RPM[i]=0.0; - else { - RPM[i] = (GMT[i] - WMD[i]) / GMT[i]; - if (RPM[i]>1.0) RPM[i]=1.0; - if (RPM[i]<0.0) RPM[i]=0.0; - } - } - } - -} - - diff --git a/cat_vol_pbtp2.m b/cat_vol_pbtp2.m deleted file mode 100644 index c0673d7f..00000000 --- a/cat_vol_pbtp2.m +++ /dev/null @@ -1,26 +0,0 @@ -%cat_vol_pbtp2 Projection-Based Thickness. -% Main thickness projection function (Dahnke et al., 2011)). -% -% [GMT,RPM] = cat_vol_pbtp2(SEG,WMD,CSFD[,opt]) -% -% GMT (3D single) .. thickness image -% RPM (3D single) .. radial position map -% SEG (3D single) .. segment image with low and high boundary bd -% (default 1=CSF, 2=GM, 3=WM) -% WMD (3D single) .. CSF distance map -% CSFD (3D single) .. CSF distance map -% -% opt .. MATLAB structure -% .bd (1x2 single) .. [low,high] boundary values (default 1.5 and 2.5) -% .CSFD .. calculate CSFD -% .PVE .. use PVE information (0=none,1=fast,2=exact) -% -% See also cat_vbdist, cat_vol_eidist, cat_vol_pbtp, compile. -% ______________________________________________________________________ -% -% Christian Gaser, Robert Dahnke -% Structural Brain Mapping Group (https://neuro-jena.github.io) -% Departments of Neurology and Psychiatry -% Jena University Hospital -% ______________________________________________________________________ -% $Id$ diff --git a/cat_vol_pbtp2.mexa64 b/cat_vol_pbtp2.mexa64 deleted file mode 100755 index 3cb6bd6f..00000000 Binary files a/cat_vol_pbtp2.mexa64 and /dev/null differ diff --git a/cat_vol_pbtp2.mexmaca64 b/cat_vol_pbtp2.mexmaca64 deleted file mode 100755 index aeb5b0c6..00000000 Binary files a/cat_vol_pbtp2.mexmaca64 and /dev/null differ diff --git a/cat_vol_pbtp2.mexmaci64 b/cat_vol_pbtp2.mexmaci64 deleted file mode 100755 index dc83d782..00000000 Binary files a/cat_vol_pbtp2.mexmaci64 and /dev/null differ diff --git a/cat_vol_pbtp2.mexw64 b/cat_vol_pbtp2.mexw64 deleted file mode 100755 index e3175433..00000000 Binary files a/cat_vol_pbtp2.mexw64 and /dev/null differ diff --git a/compile.m b/compile.m index 92e2b8a3..7cd0dcdb 100644 --- a/compile.m +++ b/compile.m @@ -85,7 +85,6 @@ 'cat_vol_simgrow.c' 'cat_vol_localstat.c' 'cat_vol_pbtp.c' - 'cat_vol_pbtp2.c' 'cat_vol_interp3f.cpp' 'cat_vol_eidist.c' 'cat_vol_genus0.c genus0.c' @@ -389,12 +388,7 @@ [d{ni},dpp] = cat_vol_pbtp(d1+1,dw,dc); r(ni) = rms(d{ni}(d1==1)) - 5.5; s(ni) = r(ni)<0.05; - % pbt2 - n{ni} = 'cat_vol_pbtp2'; - [d{ni},dpp] = cat_vol_pbtp2(d1+1,dw,dc); - r(ni) = rms(d{ni}(d1==1)) - 5.5; - s(ni) = r(ni)<0.05; - + %% test interpolation invariance % - less difference in simple structures diff --git a/mexmaca/cat_vol_pbtp2.mex b/mexmaca/cat_vol_pbtp2.mex deleted file mode 100755 index 19e6b020..00000000 Binary files a/mexmaca/cat_vol_pbtp2.mex and /dev/null differ diff --git a/mexmaci/cat_vol_pbtp2.mex b/mexmaci/cat_vol_pbtp2.mex deleted file mode 100755 index f39dc10d..00000000 Binary files a/mexmaci/cat_vol_pbtp2.mex and /dev/null differ