diff --git a/CHANGES.txt b/CHANGES.txt index bbd4bbcc..89e2c9e5 100644 --- a/CHANGES.txt +++ b/CHANGES.txt @@ -1,4 +1,30 @@ ------------------------------------------------------------------------ +r2535 | robertdahnke | 2024-02-08 00:06:27 + +Changed paths: + M CHANGES.txt +MM cat_io_xml2csv.m + M cat_main1639.m +M cat_main_updateSPM1639.m +MM cat_run.m + M cat_run_job.m + M cat_run_job1639.m +M cat_surf_createCS2.m + M cat_surf_createCS3.m +MM cat_vol_pbtsimple.m +AD cat_vol_pbtsimple_fine.m +M cat_vol_set_com.m + +1) New: Added final CSV report to cat_run with changes in cat_io_xml2csv. + Corrected bugs in cat_run and added surface values. +1) Changed: Optimized command line output of cat_vol_set_com. +3) Changed: Slight modification of the call of the surface collision correction + in cat_surf_createCS2/3. Slight changes in cat_surf_createCS2/3 for + tests of the reconstructed surfaces. +4) Changed: Various modifications in PBTsimple (see issue #19). +5) Changed: Small visual modification in cat_vol_set_com. +6) Changed: Updated alert thresholds in cat_main_updateSPM1639. +7) Changed: Tiny correction in the APP try-catch block in cat_run_job[1639].------------------------------------------------------------------------ r2534 | gaser | 2024-02-04 18:02:09 Changed paths: diff --git a/cat_io_xml2csv.m b/cat_io_xml2csv.m index ce502069..fae54079 100644 --- a/cat_io_xml2csv.m +++ b/cat_io_xml2csv.m @@ -51,6 +51,7 @@ def.avoidfields = {'help'; 'catlog'}; def.delimiter = ','; def.outdir = {''}; + def.conclusion = 1; def.dimlim = 256; % extend in case of catROI below def.report = 'default'; def.verb = 1; @@ -62,6 +63,8 @@ % - setup for CSV / TSV export? % - write error in case of too large matrices % - dependencies + % - RD20240129: handling of missing files "$FILENAME MISSED XML - FAILED PROCESSING"? + @@ -152,11 +155,14 @@ relevant_catreport_fields = .... {'qualityratings.IQR'; 'qualityratings.NCR'; 'qualityratings.ICR'; 'qualityratings.res_ECR'; 'qualityratings.res_RMS'; ... voxel-based QC measures 'software.version_cat'; 'software.revision_cat'; 'software.version_spm'; ... - ... 'qualitymeasures.SurfaceEulerNumber'; 'qualitymeasures.SurfaceDefectArea' ; ... 'subjectmeasures.vol_abs_CGW(01)'; 'subjectmeasures.vol_abs_CGW(02)'; 'subjectmeasures.vol_abs_CGW(03)'; ... 'subjectmeasures.vol_rel_CGW(01)'; 'subjectmeasures.vol_rel_CGW(02)'; 'subjectmeasures.vol_rel_CGW(03)'; ... - 'subjectmeasures.dist_thickness_kmeans'; 'subjectmeasures.surf_TSA'; 'subjectmeasures.vol_TIV'}; - + 'subjectmeasures.dist_thickness{01}(01)'; 'subjectmeasures.dist_thickness{01}(02)'; + 'subjectmeasures.surf_TSA'; 'subjectmeasures.vol_TIV'; ... + 'qualitymeasures.SurfaceEulerNumber'; 'qualitymeasures.SurfaceDefectArea'; ... + 'qualitymeasures.SurfaceIntensityRMSE'; 'qualitymeasures.SurfacePositionRMSE'; ... + 'qualitymeasures.SurfaceSelfIntersections'; ... + }; switch job.report case 'paraonly' @@ -208,12 +214,26 @@ %% extract fieldnames from structure to build a table [hdr,tab] = cat_io_struct2table(xml,fieldnames,0); + % create average + existxml = cellfun(@exist,job.files); % only for existing files + if job.conclusion + avg = cell(1,size(tab,2)); + for ci = 1:size(tab,2) % for each column + if isnumeric(cell2mat(tab(existxml>0,ci))) % for all numberic fields + avg{1,ci} = mean( cell2mat(tab(existxml>0,ci)) ); % average existing + else + avg{1,ci} = ''; + end + end + end % add index fieldnames = ['filenames'; fieldnames]; hdr = [{'filename'} hdr]; - tab = [job.files tab]; - + tab = [job.files(existxml>0) tab]; + if job.conclusion + avg = [{sprintf('average (%d of %d)',sum(existxml>0),numel(existxml))} avg]; + end % cleanup some fields ROInamelim = 30; @@ -235,12 +255,17 @@ ROI = ''; end + % header hdr{hi} = cat_io_strrep(hdr{hi},{'(','{','['},''); hdr{hi} = cat_io_strrep(hdr{hi},{')','}',']'},'_'); if hdr{hi}(end)=='_', hdr{hi}(end) = []; hdr{hi} = [hdr{hi} ROI]; end - end - table = [hdr;tab]; + end + if job.conclusion + table = [hdr;tab;avg]; + else + table = [hdr;tab]; + end @@ -267,6 +292,7 @@ end end + % write file cat_io_csv(fname,table,'','',struct('delimiter',job.delimiter,'komma','.')) diff --git a/cat_main1639.m b/cat_main1639.m index 630c0318..92290b34 100644 --- a/cat_main1639.m +++ b/cat_main1639.m @@ -799,7 +799,7 @@ [Yth1, S, Psurf, qa.createCS] = ... cat_surf_createCS3(VT,VT0,Ymix,Yl1,YMF,YT,Yb0,struct('trans',trans,'reduce_mesh',job.extopts.reduce_mesh,... required for Ypp output - 'outputpp',job.output.pp,'surf_measures',job.output.surf_measures, ... + 'outputpp',job.output.pp,'surf_measures',job.output.surf_measures, ... 'skip_registration', 1, ... 'interpV',job.extopts.pbtres,'pbtmethod',job.extopts.pbtmethod,'SRP', mod(job.extopts.SRP,10), ... 'scale_cortex', job.extopts.scale_cortex, 'add_parahipp', job.extopts.add_parahipp, 'close_parahipp', job.extopts.close_parahipp, .... 'Affine',res.Affine,'surf',{surf},'pbtlas',job.extopts.pbtlas, ... % pbtlas is the new parameter to reduce myelination effects diff --git a/cat_main_updateSPM1639.m b/cat_main_updateSPM1639.m index 7d244199..881bfd09 100644 --- a/cat_main_updateSPM1639.m +++ b/cat_main_updateSPM1639.m @@ -63,11 +63,33 @@ clear Pti; % RD20221102: Add warnings % These threshold are arbitrary choosen and updates are maybe required! + % RD202401: updated (worst cases with *) + % BUSS Median: GM=0.75, WM=0.31, CSF=0.02** + % 4397 Median: GM=0.31, WM=0.63, CSF=0.09* + % NISALS Median: GM=0.22, WM=0.25, CSF=0.01* + % ADHD Median: GM=0.84, WM=0.97, CSF=0.00* + % OASIS Median: GM=0.22, WM=0.31, CSF=0.02* + % Colins Median: GM=0.80, WM=0.97, CSF=0.06 + % BWPT Median: GM=0.96, WM=1.00, CSF=0.05 + % HR075 Median: GM=0.75, WM=0.98, CSF=0.03 + % <0.50 <0.50 <0.02 + % + % BUSS Mean: GM=0.54, WM=0.48, CSF=0.13 + % 4397 Mean: GM=0.44, WM=0.53, CSF=0.27 + % NISALS Mean: GM=0.43, WM=0.46, CSF=0.12 + % ADHD Mean: GM=0.59, WM=0.66, CSF=0.06 + % OASIS Mean: GM=0.37, WM=0.45, CSF=0.15 + % Colins Mean: GM=0.56, WM=0.72, CSF=0.29 + % BWPT Mean: GM=0.72, WM=0.83, CSF=0.37 + % HR075 Mean: GM=0.54, WM=0.79, CSF=0.20 + % <0.5 <0.6 <0.13 + if job.extopts.expertgui>1 && ... - ( any( res.spmP0.md(1:3) < [.5 .9 .3],2) || any( res.spmP0.mn(1:3) < [.4 .7 .2],2) || ... + ( any( res.spmP0.md(1:3) < [.4 .4 .02],2) || ... + any( res.spmP0.mn(1:3) < [.4 .5 .10],2) || ... any( any( res.spmP0.hstsumi(1:3,round(hbuckets*[.25,.50,.75])) < [0.6 0.5 0.4; 0.8 0.7 0.6; 0.5 0.2 0.05 ] )) ) cat_io_addwarning('cat_main_updateSPM1639:lowSPMseg',... - sprintf(['Low SPM segmentation quality with larger areas of mixed tissues. \\\\n' ... + sprintf(['Low SPM segmentation quality with larger areas of mixed tissues (smaller=worse). \\\\n' ... ' Median: GM=%0.2f, WM=%0.2f, CSF=%0.2f \\\\n' ... ' Mean: GM=%0.2f, WM=%0.2f, CSF=%0.2f'], ... res.spmP0.md(1:3),res.spmP0.mn(1:3)),2,[1 2]); @@ -77,6 +99,15 @@ % Besides the probabilty we may check also the intensities. If two % classes show a high overlap of intensities the segmentation is mabe % inaccurate. + % BUZZ CG=0.81,GW=0.38,CW=0.82,CGW=0.25 + % 4397 CG=0.72,GW=0.43,CW=0.84,CGW=0.29 + % NISALS CG=0.41,GW=0.44,CW=0.62,CGW=0.18 + % ADHD CG=0.27,GW=0.31,CW=0.41,CGW=0.17 + % OASIS CG=0.39,GW=0.37,CW=0.58,CGW=0.14 + % Colins CG=0.45,GW=0.34,CW=0.48,CGW=0.24 + % BWPT CG=0.39,GW=0.23,CW=0.42,CGW=0.16 + % HR075 CG=0.67,GW=0.50,CW=0.72,CGW=0.32 + % <0.3 <0.3 <0.3 <0.15 minprob = 32; hbuckets = 30; trange = [min(res.mn(res.lkp<6)) max(res.mn(res.lkp<6))]; trange = trange + [ -diff(trange) +diff(trange) ] / 4; @@ -98,7 +129,7 @@ res.spmP0.help = [res.spmP0.help; { sprintf('The hsti is the histogram of the T1 image for each class (prob>%0.2f,%d buckets) ',minprob/255,hbuckets)}; 'that are combined into the average intensity difference between classes histidiff with 1 for the ideal case without overlap and 0 for high overlap. ']; - if job.extopts.expertgui>1 && res.spmP0.hstidiff<0.5 + if job.extopts.expertgui>1 && res.spmP0.hstidiff<0.15 cat_io_addwarning('cat_main_updateSPM1639:highIntOverlapBetweenClasses',... sprintf('High overlap of image intensies in different classes\\\\n CG=%0.2f,GW=%0.2f,CW=%0.2f,CGW=%0.2f',res.spmP0.hstidiffCGW,res.spmP0.hstidiff),2,[1 2]); end diff --git a/cat_run.m b/cat_run.m index 5c06f9c1..089081fe 100644 --- a/cat_run.m +++ b/cat_run.m @@ -331,6 +331,7 @@ GMC = GMC ./ repmat( max(1,sum(GMC,2)) , 1 , 3); % make bright values darker color = @(QMC,m) QMC(max(1,min(size(QMC,1),round(((m-1)*3)+1))),:); colorgmt = @(GMC,m) GMC(max(1,min(size(GMC,1),round(((m-0.5)*10)+1))),:); + colorsurf = @(SI,m) SI(max(1,min(size(SI,1),round(((m-0.06)*4000)+1))),:); rps2mark = @(rps) min(10.5,max(0.5,10.5 - rps / 10)) + isnan(rps).*rps; % not used %mark2rps = @(mark) min(100,max(0,105 - mark*10)) + isnan(mark).*mark; @@ -455,11 +456,23 @@ else catgmt = {'unknown'}; end + %% surface intensity / position RMSE + cati = find(cellfun('isempty',strfind(txt,'Surface intensity / position RMSE: '))==0,1,'last'); + if ~isempty(cati) + cathd = textscan( txt{cati} ,'%s','Delimiter',':'); + cathd = textscan( cathd{1}{2},'%s','Delimiter',' '); + catSRMSE = [cathd{1}(1) cathd{1}(3)]; + else + catSRMSE = {'unknown'}; + end + + + %% search WARNINGs and ERRORs - cati = find(cellfun('isempty',strfind(txt(catis:end),'ALERT '))==0); + cati = find(cellfun('isempty',strfind(txt(catis(end):end),'ALERT '))==0); catalerts = numel(cati); - cati = find(cellfun('isempty',strfind(txt(catis:end),'WARNING '))==0); + cati = find(cellfun('isempty',strfind(txt(catis(end):end),'WARNING '))==0); catwarnings = numel(cati); @@ -598,7 +611,7 @@ cat_io_cprintf([0 0 0],idx); cat_io_cprintf([0 0 0],catlogt); if isempty(caterr) - cat_io_cprintf([0 0 0],sprintf(' % 3d.%02d minutes, ',cattime')); + cat_io_cprintf([0 0 0],sprintf('% 5d.%02d minutes, ',cattime')); % add IQR col = color(QMC,rps2mark( str2double( catiqr{1}(1:end-1) ))); @@ -636,6 +649,17 @@ if job.extopts.expertgui > 0 && ~strcmp(catgmt{1},'unknown') cat_io_cprintf(kcol,', '); cat_io_cprintf(col,sprintf('GMT=%s',strrep(catgmt{1},'%','%%'))); end + + % surf vals + if job.extopts.expertgui > 0 && ~strcmp(catgmt{1},'unknown') + colorsurf = @(SI,m) SI(max(1,min(size(SI,1),round((max(0,m-0.06)*1000)+1))),:); + cat_io_cprintf(kcol,', '); + cat_io_cprintf(colorsurf(GMC,str2double( catSRMSE{1} )),sprintf('surf=%s',catSRMSE{1})); + cat_io_cprintf(kcol,', '); + cat_io_cprintf(colorsurf(GMC,str2double( catSRMSE{2} )),sprintf('%s' ,catSRMSE{2})); + end + + % warnings allcatwarnings = allcatwarnings + catwarnings; @@ -659,11 +683,11 @@ cat_io_cprintf(kcol,', '); cat_io_cprintf(col,sprintf('%d alerts',catalerts)); end - cat_io_cprintf('n',' '); % to avoid color bug? - fprintf(' \n'); else - cat_io_cprintf(err.color,sprintf('%s\n',err.txt)); + cat_io_cprintf(err.color,sprintf('%s',err.txt)); end + cat_io_cprintf(kcol,'. '); % to avoid color bug? + fprintf(' \n'); end end cat_progress_bar('Set', cid ); @@ -672,6 +696,16 @@ end err.missed = max(0,sum( numel(job_data) ) - (err.warn0 + err.warn1 + err.warn2 + err.aff + err.vbm + err.sbm)); + + % RD202401: create CSV report for all processed files + % besides this a nice para report could be useful that include + % all fields and not just the most relevant on listed in the + % report. + try + cat_run_createCSVreport(job,BIDSfolder); + end + + %% final report fprintf('_______________________________________________________________\n'); fprintf('Conclusion of %d cases: \n',numel(job_data)); @@ -727,6 +761,7 @@ varargout{1} = jobl.vout; end + % clear useprior option to ensure that option is set back to default % for next processings cat_get_defaults('useprior',''); @@ -735,6 +770,38 @@ varargout{1} = cat_io_checkdepfiles( varargout{1} ); return %_______________________________________________________________________ +function cat_run_createCSVreport(job,BIDSfolder) +%% create a final csv with values from the XML reports + + % define input XMLs + matlabbatch{1}.spm.tools.cat.tools.xml2csv.files = job.data; + for fi = 1:numel(job.data) + [~, reportfolderfi] = cat_io_subfolders(job.data{fi}, job); + matlabbatch{1}.spm.tools.cat.tools.xml2csv.files{fi} = spm_file(job.data{fi}, ... + 'path', fullfile( spm_fileparts(job.data{fi}), reportfolderfi ), ... + 'prefix', 'cat_', 'ext', '.xml'); + end + + % RD20240129: HOW TO HANDLE MISSED FILES ? >> handle in cat_io_xml2csv + + % filename with date + [~,pp1,pp2] = spm_fileparts(BIDSfolder); + date = ''; %['_' datetime('now','Format','yyyyMMddHHmm')]; + if ~isempty(pp1), pp1 = ['_' pp1 strrep(pp2,'_','-')]; end + matlabbatch{1}.spm.tools.cat.tools.xml2csv.fname = ... + sprintf('CATxml%s%s.csv', pp1, date); + matlabbatch{1}.spm.tools.cat.tools.xml2csv.outdir = ... + {fullfile( spm_fileparts(job.data{fi}), spm_str_manip(reportfolderfi,'h'))}; + matlabbatch{1}.spm.tools.cat.tools.xml2csv.fieldnames = {' '}; + matlabbatch{1}.spm.tools.cat.tools.xml2csv.avoidfields = {''}; + matlabbatch{1}.spm.tools.cat.tools.xml2csv.report = 'default'; + + % run SPM batch + warning off; + spm_jobman('run',matlabbatch); + warning on; +return +%_______________________________________________________________________ function job = update_job(job) % set GUI specific parameter if available diff --git a/cat_run_job.m b/cat_run_job.m index 48e7158a..059d0b08 100644 --- a/cat_run_job.m +++ b/cat_run_job.m @@ -566,8 +566,8 @@ function cat_run_job(job,tpm,subj) if job.extopts.APP == 1070 && ~ppe.affreg.highBG && ... ( ~isfield(job,'useprior') || isempty(job.useprior) ) stime = cat_io_cmd('Affine preprocessing (APP)'); + Ysrc = single(obj.image.private.dat(:,:,:)); try - Ysrc = single(obj.image.private.dat(:,:,:)); [Ym,Yt,Ybg,WMth] = cat_run_job_APP_init1070(Ysrc,vx_vol,job.extopts.verb); %#ok catch apperr %% very simple affine preprocessing ... only simple warning diff --git a/cat_run_job1639.m b/cat_run_job1639.m index 74a589b8..29fe9de8 100644 --- a/cat_run_job1639.m +++ b/cat_run_job1639.m @@ -541,8 +541,8 @@ function cat_run_job1639(job,tpm,subj) if job.extopts.APP == 1070 && ~ppe.affreg.highBG && ... ( ~isfield(job,'useprior') || isempty(job.useprior) ) stime = cat_io_cmd('Affine preprocessing (APP)'); + Ysrc = single(obj.image.private.dat(:,:,:)); try - Ysrc = single(obj.image.private.dat(:,:,:)); [Ym,Yt,Ybg,WMth] = cat_run_job_APP_init1070(Ysrc,vx_vol,job.extopts.verb); %#ok catch apperr %% very simple affine preprocessing ... only simple warning diff --git a/cat_surf_createCS2.m b/cat_surf_createCS2.m index 04f3a515..0919a109 100644 --- a/cat_surf_createCS2.m +++ b/cat_surf_createCS2.m @@ -1333,14 +1333,15 @@ %% call collision correction % RD202108: Use further iterations if self-intersections are still very high. % (test data was an high resolution ex-vivo chimp PD image that had still strong SIs after first correction) - SIOs = 100; SIs = 80; maxiter = 2; iter = 0; - while SIs>5 && SIs5 && SIs=2,'verb',verblc,'mat',Smat.matlabIBB_mm,'vx_vol',vx_vol)); + [CS,facevertexcdata,SIs] = cat_surf_fun('collisionCorrectionPBT',CS,facevertexcdata,Ymfs,Yppi,... + struct('optimize',iter<2 && opt.SRP>=2,'verb',verblc,'mat',Smat.matlabIBB_mm,'vx_vol',vx_vol)); if verblc, fprintf('\b\b'); end if strcmpi(spm_check_version,'octave') && iter == 1 cat_io_addwarning('cat_surf_createCS2:nofullSRP','Fine correction of surface collisions is not yet available under Octave.',2) - else + elseif iter == 1 % to keep it fast we just do this once [CS,facevertexcdata,SIs] = cat_surf_fun('collisionCorrectionRY' ,CS,facevertexcdata,Ymfs,struct('Pcs',Pcentral,'verb',verblc,'mat',Smat.matlabIBB_mm,'accuracy',1/2^3)); end end @@ -1512,7 +1513,8 @@ cat_io_FreeSurfer('write_surf_data',Pthick,facevertexcdata); % final surface evaluation - res.(opt.surf{si}).createCS_final = cat_surf_fun('evalCS',loadSurf(Pcentral),cat_io_FreeSurfer('read_surf_data',Ppbt),facevertexcdata,Ymfs,Yppi,Pcentral,Smat.matlabIBB_mm,debug,cat_get_defaults('extopts.expertgui')>1); +cat_surf_fun('saveico',CS,facevertexcdata,Pcentral,sprintf('createCS_3_collcorr_%0.2fmm_vdist%0.2fmm',opt.interpV,opt.vdist),Ymfs,Smat.matlabIBB_mm); + res.(opt.surf{si}).createCS_final = cat_surf_fun('evalCS',loadSurf(Pcentral),cat_io_FreeSurfer('read_surf_data',Ppbt),facevertexcdata,Ymfs,Yppi,Pcentral,Smat.matlabIBB_mm,debug + (cat_get_defaults('extopts.expertgui')>1),cat_get_defaults('extopts.expertgui')>1); else % otherwise simply copy ?h.pbt.* to ?h.thickness.* copyfile(Ppbt,Pthick,'f'); diff --git a/cat_surf_createCS3.m b/cat_surf_createCS3.m index 9491f076..8a86c197 100644 --- a/cat_surf_createCS3.m +++ b/cat_surf_createCS3.m @@ -78,7 +78,8 @@ def.outputpp.native = 0; % output of Ypp map for cortical orientation in EEG/MEG def.outputpp.warped = 0; def.outputpp.dartel = 0; - + def.skip_registration = 0; + % options that rely on other options opt = cat_io_updateStruct(def,opt); clear def; opt.vol = any(~cellfun('isempty',strfind(opt.surf,'v'))); % only volume-based thickness estimation @@ -86,7 +87,7 @@ opt.interpV = max(0.1,min([opt.interpV,1.5])); % general limitation of the PBT resolution % for testing only - skip_registration = debug & 1; + skip_registration = debug | opt.skip_registration; % we should test this more detailed! % opt.add_parahipp = 0; @@ -691,9 +692,16 @@ %cmd = sprintf('CAT_MarchingCubesGenus0 -fwhm "3" -thresh "%g" "%s" "%s"',th_initial,Vpp_side.fname,Pcentral); % -pre-fwhm = -1: masked smoothing with 1mm % -post-fwhm = 2: psot smoothing with 2mm (use lower threshold of 0.49) - cmd = sprintf('CAT_VolMarchingCubes -pre-fwhm "-1" -post-fwhm "2" -thresh "%g" "%s" "%s"',th_initial,Vpp_side.fname,Pcentral); - cat_system(cmd,opt.verb-3); - + cmd = sprintf('CAT_VolMarchingCubes -pre-fwhm "-1" -post-fwhm "1" -thresh "%g" "%s" "%s"',th_initial,Vpp_side.fname,Pcentral); + cat_system(cmd,opt.verb-3); + + % Collins-without: 2.5996 ± 0.6292 mm, 0.0798 / 0.1096, 10.29% (121.38 cm²) 24.19% (285.46 cm²) + % Collins-with: 2.5713 ± 0.6525 mm, 0.0723 / 0.0934, 8.51% (98.93 cm²) 23.79% (276.42 cm²) + cmd = sprintf(['CAT_DeformSurf "%s" none 0 0 0 "%s" "%s" none 0 1 -1 .1 ' ... + 'avg -0.1 0.1 .2 .1 5 0 "0.5" "0.5" n 0 0 0 %d %g 0.0 0'], ... + Vpp_side.fname,Pcentral,Pcentral,50,0.001); + cat_system(cmd,opt.verb-3); + end % @@ -710,15 +718,6 @@ % evaluate and save results if isempty(stime), stime = clock; end fprintf('%5.0fs',etime(clock,stime)); stime = []; if 1, fprintf('\n'); end %debug - %{ - res.(opt.surf{si}).createCS_init = cat_surf_fun('evalCS',CS,cat_surf_fun('isocolors',CS,Yth1i,Smat.matlabIBB_mm),[],Ymfs,Yppi,Pcentral,Smat.matlabIBB_mm,debug); - if 0 %debug - % save surface for further evaluation - cat_surf_fun('saveico',CS,cat_surf_fun('isocolors',Yth1i,CS.vertices,Smat.matlabIBB_mm),Pcentral,sprintf('createCS_1_init_pbtres%0.2fmm',opt.interpV),Ymfs,Smat.matlabIBB_mm); - else - fprintf('\n'); - end - %} end % for use of average prior in long. pipeline we have to deform the average mesh to current pp distance map @@ -779,7 +778,7 @@ stime = cat_io_cmd(' Reduction of surface collisions with optimization:','g5','',opt.verb,stime); end CS = loadSurf(Pcentral); - SIOs = 100; SIs = 80; maxiter = 1; iter = 0; verblc = debug; % maxiter =2 + SIOs = 100; SIs = 80; maxiter = 2; iter = 0; verblc = debug; while SIs>5 && SIs1, fprintf('\n'); end - if debug + if debug || cat_get_defaults('extopts.expertgui')>1 cat_surf_fun('saveico',CS,cat_surf_fun('isocolors',Yth1i,CS.vertices,Smat.matlabIBB_mm),Pcentral,'',Ymfs,Smat.matlabIBB_mm); end - res.(opt.surf{si}).createCS_final = cat_surf_fun('evalCS',loadSurf(Pcentral),cat_io_FreeSurfer('read_surf_data',Ppbt),facevertexcdata,Ymfs,Yppi,Pcentral,Smat.matlabIBB_mm,debug,cat_get_defaults('extopts.expertgui')>1); + res.(opt.surf{si}).createCS_final = cat_surf_fun('evalCS',loadSurf(Pcentral),cat_io_FreeSurfer('read_surf_data',Ppbt),facevertexcdata,Ymfs,Yppi,Pcentral,Smat.matlabIBB_mm,debug + (cat_get_defaults('extopts.expertgui')>1),cat_get_defaults('extopts.expertgui')>1); else % otherwise simply copy ?h.pbt.* to ?h.thickness.* copyfile(Ppbt,Pthick,'f'); @@ -1000,7 +999,7 @@ Porthcolor = [ Porthcolor(1:end-1) '}']; Porthnames = [ Porthnames(1:end-1) '}']; - if debug + if 1 %debug fprintf(' Show surfaces in orthview: %s\n',spm_file(Po ,'link',... sprintf('cat_surf_fun(''show_orthview'',%s,''%s'',%s,%s)',Porthfiles,Po,Porthcolor,Porthnames))) ; end diff --git a/cat_vol_pbtsimple.m b/cat_vol_pbtsimple.m index 87d2a35d..81f1b783 100644 --- a/cat_vol_pbtsimple.m +++ b/cat_vol_pbtsimple.m @@ -1,24 +1,81 @@ function [Ygmt,Ypp] = cat_vol_pbtsimple(Yp0,vx_vol,opt) %cat_vol_pbtsimple. Simple cortical thickness/position estimation. -% Voxel-based distance estimation and projection-based thickness (PBT) -% and distance-based surface position estimation. Uses a label map as -% input. Required a isotropic input map. +% Estimation of voxel-based distances and projection-based thickness (PBT) +% and surface position based on a label map. Required isotropic input. +% +% After further optimisation the function becomes more complex again. +% However, the supersimple parameter can be used to switch off extra steps. +% Although, the values look quite similar difference in the surface are +% good detectable for small structures, in particular occipital sulci. % % [Ygmt, Ypp] = cat_vol_pbtsimple( Yp0 , vx_vol, opt ) % % Ygmt .. GM thickness map % Ypp .. percentage position map % Yp0 .. tissue label map (1-CSF, 2-GM, 3-WM) -% vx_vol .. voxel-size (in mm, default = 1) -% opt .. parameter structure -% .levels .. +% vx_vol .. voxel-size (in mm, default=1) + +% opt .. parameter structure % -% See also cat_vol_pbt, cat_vol_createCS3. -% ______________________________________________________________________ +% .supersimple (0-no, 1-yes; default=1) +% WARNING: Activation will run only the basic routines +% without enhanced refinements! +% Although surface intensity and position values are not bad +% the number of self-intersections indicate further problems. +% In particular, small sulci are often affected by blurring. +% +% .levels (integer; default=8) +% Number of dual distance estimations to reduce sampling effects. +% With logarithmic improvement and good results between 2 and 16. +% +% .extendedrange (0-no, 1-yes; default=1) +% Uses also the PVE range to estimate the distance there. +% Important to avoid thickness underestimations. +% +% .correctoffeset (0-none, 1-fixed, 2-adaptive; default=2) +% Correction for the offset of boundaries beyond the paired concept. +% Only minor effects. +% +% .distcleanup (0-no, 1-yes; default=1); +% Remove distance outliers, e.g., blood vessels, by the histogram but +% more important by a first very smooth thickness approximation that +% define local outliers. Although, this gives only tiny global +% improvements, the local differences are essential, e.g. occipital +% close to the cerebellum and major veins. +% +% .gyrusrecon (0-none, 1-simple, 2-complex; default=3) +% Additional calls of PBT to reconstruct also the gyri to reduce +% thickness overestimations. +% 0 - none +% Less accurate thickness with gyral overestimation but more +% accurate outer surface compared to method 1. +% 1 - simple approach +% Independent reconstruction of sulci and gyri that simply takes +% the minimum thickness. This causes light underestimations of +% the outer surface position in gyral regions (i.e. running too +% much inside) but better thickness values. +% 2 - complex approach (default) +% Stepwise reconstruction of sulci and gyri (i.e. dependent). +% This is sensitive for (blood vessel) artefacts that is handled +% by distcleanup step. +% +% .range (real value <=0.5, good between 0.2 and 0.4; default=0.3) +% Limitation of the offset of multiple thickness levels to avoid +% running into partial volume effects with thickness overestimation. % -% Dahnke, R; Yotter R; Gaser C. -% Cortical thickness and central surface estimation. -% NeuroImage 65 (2013) 226-248. +% .keepdetails (0-off, 1-sulci, 2-sulci+gyri; default=1) +% Enhance thin sulci (and gyri) to avoid blurring. +% Small global differences but important in occipial regions. +% Gyri enhancement is not optimal yet. +% +% .sharpening (0-no, 1-yes; default=1) +% Further optimization of the maps that can help to save small sulci. +% Sharpening the Ypp map to support more better resampling to lower +% resolution for the 0.5 layer adopted for slight changes that improves +% especially the position RMSE value by ~0.01 (ie, unclear if this is +% really better). +% +% See also cat_vol_pbt, cat_vol_createCS[123]. % ______________________________________________________________________ % % Christian Gaser, Robert Dahnke @@ -27,363 +84,489 @@ % Jena University Hospital % ______________________________________________________________________ % $Id$ - -%#ok<*UNRCH> + +%#ok<*UNRCH,*HIST> if ~exist('opt','var'), opt = struct(); end + + if 0 + % test old function that resulted in better intensity and position values + [Ygmt,Ypp] = cat_vol_pbtsimple_goodold(Yp0,vx_vol,opt); + return + end - def.WMTC = 1; % Partial voxel-based topology correction of small - % WM-defects based on the volume output of cat_vol_genus0 - def.levels = 4; % Number of dual distance estimate (requires refinement) with 0 - % for the classic approach with 1.5 and 2.5 as CSF and WM boudnary. - % Good values are between 2 and 4, whease higher values (>6) - % run into problems for interpolation artifacts of close to the - % full tissue class. - def.refine = 1; % Refined WM distance based on CSF distance (myelination correction) - def.gyrusrecon = 1; % additional PBT gyri reconstruction (TRUE/FALSE) - % this reduce thickness overestimations but maybe underestimates - % the outer surface position in gyral regions (running to much inside) - def.correctoffeset = 2; % not really required if no refinement is used - def.extendedrange = 1; % may causes closed gyri - def.sharpening = 0; % sharpening the Ypp map to support more better resampling to lower resolution for the 0.5 layer - worse - def.keepdetails = 1; % enhance thin sulci to avoid blurring (0-off, 1-sulci, 2-sulci+gyri) - opt = cat_io_checkinopt(opt,def); + % Default parameter settings with short impression for some Collins + % processing to roughly quantify the global effect. + % + % Values: + % thickness (mn±sd), surface intensity / position RMSE, SI=self-intersections + % + % Final tests with the following subjects are important : + % Collins, HR075, 4397, OASIS031 (WMHs), BUSS02 (Child), + % BWPT (Phantom), NISALS_UTR (PVEs), ADHD200NYC (LowQuality/Motion) + % + % Overall, it is important to eep an eye on + % * thickness values >> BWPT + % * unterestimation of very thick young cortices >> BUSS02 + % * local breaks/defects with thickness underestimation in older subjects with WMHs >> OASIS31, NISALS_UTR + % * or in case of motion artefact >> ADHD200NYC + % + def.supersimple = 0; % Refined WM distance based on CSF distance (myelination correction) + % This internal options is to turn of all feature functions. + % Collins-1: 2.5713 ± 0.6525 mm, 0.0723 / 0.0934, SI=16.15% + % Collins-0: 2.5633 ± 0.6876 mm, 0.0728 / 0.0805, SI= 7.67% + % -0.0005 / 0.0129, +8.48% + % + def.WMTC = 0; % Voxel-based topology correction of small WM-defects + % based on the volume output of cat_vol_genus0 + % ... not really working ... just keep it now for sentimental reasons before it will be remove in future + % + def.levels = 8; % Number of dual distance estimates. + % Larger values are more accurate but need more time (log-change). + % Good values are between 2 and 16 with limited range parameter. + % + def.correctoffeset = 2; % correct for layer intensity (0-no, 1-yes-simple, 2-yes-complex; tiny improvement) + % Collins-0: 2.5647 ± 0.6883 mm, 0.0726 / 0.0814, SI=7.78% + % Collins-1: 2.5633 ± 0.6876 mm, 0.0728 / 0.0805, SI=7.67% + % -0.0002 / 0.0009, +0.11% + % + def.extendedrange = 1; % may causes closed gyri (0-no, 1-yes) + def.distcleanup = 1; % remove distance outliers, e.g., blood vessels + % (0-no, 1-yes; tiny global improvement but locally pretty important) + % Collins-1: 2.6693 ± 0.6199 mm, 0.0734 / 0.0746, SI=4.975% + % Collins-0: 2.6493 ± 0.6147 mm, 0.0728 / 0.0730, SI=4.530% + % +0.0006 /+0.0016, +0.445% + % + def.gyrusrecon = 2; % Additional PBT gyri reconstruction + % this reduce thickness overestimations but maybe underestimates + % the outer surface position in gyral regions (running too much inside) + % 0 - none 0.0666 / 0.0856, SI=10.19% + % 1 - simple old approach 0.0677 / 0.0838, SI=17.09% + % 2 - more complex approach 0.0729 / 0.0745, SI= 4.30% *** + % + def.range = 0.3; % Default value for range extension (should be between 0.2-0.4) + % + def.keepdetails = 1; % enhance thin (occipital) sulci (and gyri) to avoid blurring + % (0-no, 1-yes; 2-yes(+gyri) worse values but pretty important! + % (0-off, 1-sulci, 2-sulci+gyri) + % Collins-0: 2.6413 ± 0.6082 mm, 0.0708 / 0.0722, SI=4.585% + % Collins-1: 2.6509 ± 0.6151 mm, 0.0728 / 0.0730, SI=5.365% + % -0.0020 /-0.0012, -1.220% + % + def.sharpening = 1; % sharpening the Ypp map to avoid blurring while + % resampling to lower resolution + % (0-no, 1-yes; reduce collisions) + % Collins-0: 2.5654 ± 0.6841 mm, 0.0730 / 0.0845, SI=8.22% + % Collins-1: 2.5633 ± 0.6876 mm, 0.0728 / 0.0805, SI=7.67% + % 0.0002 / 0.0040, SI=0.55% + % + opt = cat_io_checkinopt(opt,def); -% == preparations == -% should maybe be better part of create CS + % just a fast option to switch of the extra functions + if opt.supersimple % avoid extras + opt.WMTC = 0; % WM topology correction + opt.levels = 0; % number of dual distance measurements + opt.correctoffset = 0; % correct of dual displacements (in principle not required) + opt.extendedrange = 0; % use further PVE area + opt.distcleanup = 0; % important correction of local outliers like blood vessels + opt.gyrusrecon = 0; % reconstruct gyri (in very young/old subjects with thin WM) + opt.keepdetails = 0; % important for small (occipital) sulci + opt.sharpening = 0; % reduce blurring of the Ypp map + end - if opt.WMTC - % RD20231224: for WM topology smoothing - % the idea is to apply the voxel-based correction and _close_ small WM holes - % (<10 voxel) for the 2.75 and 2.25 boundaries for values above 2. + % display paras for developers + if cat_get_defaults('extopts.expertgui') > 1 + fprintf('\nPBTsimple parameters:\n'); + disp(opt); + end - % indicate isolated holes and replace by median of the neighbors - Yp0(Yp0<0.35 & ~cat_vol_morph(Yp0<1,'l')) = 1; % close major wholes in the WM - Ymsk = Yp0==-1 & cat_vol_morph(Yp0>0.9,'d',1); % filter small wholes close to the WM - Yp0 = cat_vol_median3(single(Yp0)+1,Ymsk,~Ymsk)-1; - % indicate isolated objects and replace by median of the neighbors - Yp0(Yp0>0.65 & cat_vol_morph(Yp0==-1,'l')) = -1; - Ymsk = Yp0>0.95 & cat_vol_morph(Yp0<-0.9,'d',1); - Yp0 = cat_vol_median3(single(Yp0)+1,Ymsk,~Ymsk)-1; - % RD20210401: but there are some new background dots - Ymsk = Yp0==-1 & cat_vol_morph(Yp0>-1,'lc'); % close major wholes in the WM - Yp0 = cat_vol_median3(single(Yp0)+1,Ymsk,~Ymsk)-1; - end - - - - if opt.WMTC - % topology correction for WM based on the volume output of cat_vol_genus0 - - % Closing of the WM object to avoid smaller holes. - % Opening is here not useful as the WM ! - tclevels = [2.75, 2.25]; - for ii = 1:1 - for tci = 1:numel(tclevels) - evalc(sprintf('Yppc = cat_vol_genus0(Yp0, %0.2f,0);',tclevels(tci))); - %if tclevels(tci) > 2.5 - Yholes = Yppc==1 & Yp0=tclevels(tci); - Ybridges = Ybridges & ~cat_vol_morph(Ybridges,'l',[1e4,8]); % allow only small corrections - Yp0(Ybridges) = tclevels(tci) - 0.01; % open - clear Ybridges; - %end - end - end - - - % object correction in the background (opening of GM-WM objects) - tclevels = [1.75, 1.25, 1.01]; - for tci = 1:numel(tclevels) - Yp0(Yp0>tclevels(tci) & ~cat_vol_morph(Yp0>tclevels(tci), ... - 'ldo',2.5 - tclevels(tci))) = tclevels(tci) - 0.01; - end + % WM topology correction (relevant in thin WM subjects such as OASIS31): + if opt.WMTC + Yp0 = WMTC(Yp0); end + % (1) Distance estimation: % -------------------------------------------------------------------- - if opt.levels == 0 + if opt.levels == 0 % Classic CSF and WM distance maps based on a binary boundary without - % partial volume effect (only for tests). - Ycd = cat_vbdist(single(Yp0 < 1.5), Yp0 < 3); - Ywd = cat_vbdist(single(Yp0 > 2.5), Yp0 > 1); - - elseif opt.levels > 0 && opt.refine == 0 - % Simple CSF and WM distance maps with partial volume effect by multiple - % distance estimations but without further corrections. - Ycd = cat_vol_eudist( max(0,min(1,2.0 - Yp0)), Yp0 <= 2.50, opt.levels, opt.correctoffeset); - Ywd = cat_vol_eudist( max(0,min(1,Yp0 - 2.0)), Yp0 >= 1.50, opt.levels, opt.correctoffeset); - opt.extendedrange = 0; - - elseif opt.levels > 0 && opt.refine > 0 - % Enhanced CSF and WM distance maps with partial volume effect by - % multiple distance estimations but without further corrections. - [Ycd, Ywd] = cat_vol_cwdist(Yp0,opt); - + % partial volume effect (only as simplest approach and for tests). + Ycd = cat_vbdist( single(Yp0 < 1.5), Yp0 < 2.5 + opt.range*opt.extendedrange ); + Ywd = cat_vbdist( single(Yp0 > 2.5), Yp0 > 1.5 - opt.range*opt.extendedrange ); + + else + if opt.supersimple + % Simple CSF and WM distance maps with partial volume effect by multiple + % distance estimations but without further corrections. + Ycd = cat_vol_eudist( max(0,min(1,2.0 - Yp0)), ... + Yp0 <= 2.5 + opt.range*opt.extendedrange, ... + opt.levels, opt.correctoffeset); + Ywd = cat_vol_eudist( max(0,min(1,Yp0 - 2.0)), ... + Yp0 >= 1.5 - opt.range*opt.extendedrange, ... + opt.levels, opt.correctoffeset); + + else + % Enhanced CSF and WM distance maps with partial volume effect by + % multiple distance estimations but without further corrections. + [Ycd, Ywd] = cat_vol_cwdist(Yp0, vx_vol, opt); + + end end - - + % (2) Thickness estimation: % -------------------------------------------------------------------- % Typically, only the sulci are reconstructed but also thins gyral - % structures suffer from blurring by low resolution or artificats. - % The idea is to reconstruct both sulci as well as gyris and just use + % structures suffer from blurring by low resolution or artefacts. + % The idea is to reconstruct both sulci as well as gyri and just use % the minimum thickness. The missing update of the CSF distance is not % optimal. % However, there is some small evidence (ie. better intensity/position % RMSE) and surface look (i.e. less errors in topology correction) that % this still support some small benefits. % - Basic tests in ADHD200NYC and Collins. - if ~opt.gyrusrecon - % Using the PBT apporach only to reconstruct the sulci. - % As this is the more simple apprach we should keep and test it in - % case of pipeline changes. + if opt.gyrusrecon == 0 + % Using the PBT apporach only to reconstruct the sulci. + % remove highly distant outliers + [Yp0, Ywd, Ycd] = cleanupVessels(Yp0, Ywd, Ycd, opt.distcleanup); + % projection-based thickness mapping Ygmt = cat_vol_pbtp( round(Yp0) , Ywd, Ycd); + + % now correct also this values (it is a bit better this way) + Ycd = max(0,Ycd - 0.5); Ywd = max(0,Ywd - 0.5); - if opt.levels >= 0 - Ycd = max(0,Ycd - 0.5); Ywd = max(0,Ywd - 0.5); % now correct also this values - end - % minimum to reduce issues with meninges - Ygmt = min(Ygmt, Ycd + Ywd); + Ygmt = min( max(0,Ygmt - 0.5) , Ycd + Ywd); % this should avoid sulcal overestimations a bit + %Ygmt = min(Ygmt, Ycd + Ywd); - else - % Using PBT to reconstruct the sulci and the gyri. + elseif opt.gyrusrecon == 1 % EXTRA + % Using PBT to reconstruct the sulci and the gyri to get the minimal + % thickness. Classic quite simple approach. + % remove highly distant outliers + [Yp0, Ywd, Ycd] = cleanupVessels(Yp0, Ywd, Ycd, opt.distcleanup); + % reconstruct sulci as well as gyri Ygmt1 = cat_vol_pbtp(round(Yp0) , Ywd, Ycd); Ygmt2 = cat_vol_pbtp(round(4-Yp0), Ycd, Ywd); Ygmt2 = max(Ygmt2 , 1.5 .* (Ygmt2>0)); %only in thick regions - if opt.levels >= 0 - Ycd = max(0,Ycd - 0.5); Ywd = max(0,Ywd - 0.5); % now correct also this values - end - + % now correct also this values + Ycd = max(0,Ycd - 0.5); Ywd = max(0,Ywd - 0.5); + % avoid meninges ! Ygmt1 = min(Ygmt1, Ycd + Ywd); Ygmt2 = min(Ygmt2, Ycd + Ywd); % average GM thickness maps - Ygmt = min(cat(4, Ygmt1, Ygmt2),[],4); + Ygmt = min(cat(4, Ygmt1, Ygmt2),[],4); clear Ygmt1 Ygmt2 - clear Ygmt1 Ygmt2 - end + elseif opt.gyrusrecon == 2 % EXTRA + % Using PBT to reconstruct the sulci and then the gyri. + + % remove highly distant outliers + [Yp0c, Ywdc, Ycdc] = cleanupVessels(Yp0, Ywd, Ycd, opt.distcleanup); + + % 1. sulcus-reconstruction: + % first, we just use the classic estimation to get a corrected CSF distance + Ygmt1 = cat_vol_pbtp( round(Yp0c) , Ywdc, Ycdc); + Ygmt1 = cat_vol_approx(Ygmt1,'nh'); + YM = Yp0 > 1.5 - opt.range * opt.extendedrange & ... + Yp0 < 2.5 + opt.range * opt.extendedrange & ... + Ygmt1 > eps; + Ycdc = Ycd; Ycdc(YM) = min(Ycd(YM), Ygmt1(YM) - Ywd(YM)); + [~,I] = cat_vbdist(single(Yp0<2.5),Yp0<2.5 + opt.range * opt.extendedrange); + Ycdc = Ycdc(I); clear I; + Ycdc = min(Ycd, max(eps*(Ycd>0), cat_vol_median3(Ycdc, YM, true(size(Yp0)), 0.2)) ); % remove outliers + + % 2. gyrus reconstruction: + % now, we can process the inverse/dual case go get + Ygmt2 = cat_vol_pbtp(round(4 - Yp0), Ycdc, Ywd); + Ygmt2 = cat_vol_approx(Ygmt2,'nh'); + YM2 = YM & Ycdc > median(Ygmt2(:)); + Ywdc = Ywd; Ywdc(YM2) = min(Ywd(YM2), Ygmt2(YM2) - Ycdc(YM2)); clear Ygmt2; + [~,I] = cat_vbdist(single(Yp0>1.5),Yp0>1.5 - opt.range * opt.extendedrange); + Ywdc = Ywdc(I); clear I; + Ywdc = min(Ywd, max(eps*(Ywd>0), cat_vol_median3(Ywdc, YM, true(size(Yp0)), 0.2) ) ); % remove outliers + + % 3. final estimation: + % having now the correct distance values the finale thickness estimations are applied + Ygmt = cat_vol_pbtp( round(Yp0) , Ywdc, Ycdc); + Ymsk = cat_vol_smooth3X( cat_vol_morph( Yp0==1 | Ywdc > median(Ygmt(:)) + 4*std(Ygmt(:)) + 1 ,'c'),1); + Ywdc2 = Ywdc; Ywdc2(Ymsk>0.5) = 0; + Ycdc2 = Ycdc; Ycdc2(Ymsk>0.5) = 0; + Yp0c = min( max(1,3 - 2*Ymsk), Yp0); + Ygmt = cat_vol_pbtp(round(Yp0c) , Ywdc2, Ycdc2); + + % remove outliers + Ygmt( Ygmt > Ygmt1 .* (median(Ygmt1(Ygmt1(:)>0)) / median(Ygmt(Ygmt(:)>0))) ) = 0; clear Ygmt1; + + Ycd = Ycdc; Ywd = Ywdc; + end + + % addition correction based on the thickness phantom (more required in the CS3 pipeline but why) + % this could come from the interpolation + if 0 %opt.levels == 0 && vx_vol < 1 + Ycd = max(0,Ycd - 0.25); + Ywd = max(0,Ywd - 0.25); + Ygmt = max(0,Ygmt - 0.50); + end % (3) Approximation of non GM voxels for resampling: % -------------------------------------------------------------------- Ygmt = cat_vol_approx(Ygmt,'nh'); - % code for test and backup - Ygmto=Ygmt+0; Ywdo=Ywd+0; Ycdo=Ycd+0; - %% - Ygmt=Ygmto+0; Ywd=Ywdo+0; Ycd=Ycdo+0; - - % Although distances and thickness are quite good, PBT slightly trend to - % over-estimate the thickness in sulcal regions without CSF as the middle - % voxel is counted for both sides (simplified). Moreover, initial surface - % are partially created just on the original internal resolution (~1 mm), - % what result in blurring of small sucli. To keep these structures, we - % optimize regions where blurring/closing is happening by making them a - % bit thinner and correcting the CSF distance (keepdetails>0). The dual - % operation (of avoiding of blurring small gyri) can also be applied but - % is expected to have less effects as these structures are already quite - % save by the WM (keepdetails>1). - if opt.keepdetails - % estimate current percentage map (same as bellow) to identify - % and correct problematic areas - YM = Yp0 > 1.5 - 0.45 * opt.extendedrange & ... - Yp0 < 2.5 + 0.45 * opt.extendedrange & ... - Ygmt > eps; - Ycdc = Ycd; Ycdc(YM) = min(Ycd(YM), Ygmt(YM) - Ywd(YM)); - Ypp = zeros(size(Yp0),'single'); Ypp(Yp0 >= 2.5 + 0.45 * opt.extendedrange) = 1; - Ypp(YM) = Ycdc(YM) ./ (Ygmt(YM) + eps); - - % reduce thickness and CSF distance if this results in smoothing - % do this only in thin regions - for ppth = 1:-.1:0.1 - Yblurredsulci = Ypp=ppth,'c',1) & Ygmt0) = max(0,min(Ycd(Yblurredsulci>0), Ygmt(Yblurredsulci>0) - Ywd(Yblurredsulci>0))); - end - - % reduce thickness and WM distance if this results in smoothing - % ... slow, worse values, - if opt.keepdetails == 2 - for ppth = 0.1:.1:0.9 - Yblurredgyri = Ypp>=ppth & cat_vol_morph(Yppmedian(Ygmt(:))/2; - Yblurredgyri = cat_vol_smooth3X(Yblurredgyri,2); - Ygmt = max(min(0,Ygmto + 1/mean(vx_vol)), Ygmt .* min(1.5,1 + .1 .* ppth*Yblurredgyri)); - Ywd(Yblurredgyri>0) = max(0,min(Ywd(Yblurredgyri>0), Ygmt(Yblurredgyri>0) + Ycd(Yblurredgyri>0))); - end - end - % final smoothing - spm_smooth(Ygmt,Ygmt,.25/mean(vx_vol)); + + % (*) EXTRA: optimize map - enhancement of thin areas + if opt.keepdetails + [Ywd, Ycd, Ygmt] = keepdetails(Yp0, Ywd, Ycd, Ygmt, vx_vol, ... + opt.range * opt.extendedrange, opt.keepdetails); end + + % (4) Estimate percentage position map: % -------------------------------------------------------------------- % We first create a corrected CSF distance map with reconstructed sulci. % If gyri were reconstructed too than also the WMD would have to be % corrected to avoid underestimation of the position map with surfaces % running to close to the WM. - overrange = .0; % range extention of the Ypp was resulting in worse - % T1 and position values and more topology defects - YM = Yp0 > 1.5 - 0.45 * opt.extendedrange & ... - Yp0 < 2.5 + 0.45 * opt.extendedrange & ... - Ygmt > eps; - Ycdc = Ycd; Ycdc(YM) = min(Ycd(YM), Ygmt(YM) - Ywd(YM)); - Ypp = zeros(size(Yp0),'single'); Ypp(Yp0 >= 2.5+0.45*opt.extendedrange) = 1; - Ypp(YM) = Ycdc(YM) ./ (Ygmt(YM) + eps); - if overrange>0 - Yppl = -2 + 2*smooth3(Ypp>0); Yppl(Ypp>0) = 0; - Ypph = 2*smooth3(Ypp>=1); Ypph(Ypp<1) = 0; - Ypp = min(1 + overrange,max(0 - overrange, Ypp + Yppl + Ypph )); - end + YM = Yp0 > 1.5 - opt.range * opt.extendedrange & ... + Yp0 < 2.5 + opt.range * opt.extendedrange & ... + Ygmt > eps; + Ycdc = Ycd; + Ycdc(YM) = min(Ycd(YM), Ygmt(YM) - Ywd(YM)); + Ypp = zeros(size(Yp0),'single'); + Ypp(Yp0 >= 2.5 + opt.range * opt.extendedrange) = 1; + Ypp(YM) = max(0,Ycdc(YM) ./ (Ygmt(YM) + eps)); + Ypp(Ypp<.9 & Yp0>2.5) = 1; clear Ycdc YM; -%% + + % (5) Voxel-size resolution correction: % -------------------------------------------------------------------- Ygmt = Ygmt * mean(vx_vol); + - - if 0 - % FINAL TOPOLOGY CORRECTION - % RD202312: experimental idea - remove in future - % final voxel-based topology correction of some levels - % however, the topology correction in the surface create should be more powerfull - tclevels = 1:-0.25:0; - for tci = 1:numel(tclevels) - evalc(sprintf('Yppc = cat_vol_genus0(Ypp, %0.2f,0);',tclevels(tci))); - if tclevels(tci) > 0.5 - Ypp(Yppc==1 & Ypp=tclevels(tci)) = tclevels(tci) - 0.01; % open - end - end - end - - - % sharpening - % RD202312: experimental idea - remove in future - % - helpful to compensate the downsampling - % - looks nice in sulci and gyri but the values (intensity and position) - % but also the results in ADNI are worse - % - as the position slighly changes, also the thickness is effected later - % by correcting the - if opt.sharpening - Ypps = Ypp + 1*(Ypp - cat_vol_smooth3X(Ypp,1/mean(vx_vol))) + ... - 2*(Ypp - cat_vol_smooth3X(Ypp,2/mean(vx_vol))) + ... - 4*(Ypp - cat_vol_smooth3X(Ypp,4/mean(vx_vol))); - - % final smoothing to prepare surface reconstruction that also correct for WM topology issues - spm_smooth(Ypps,Ypps,.5/mean(vx_vol)); - %Ypp = min(1 + overrange,max(0 - overrange,Ypps)); % old version - % New version that only change mid-position values relevant for surface - % creation. Although this avoid the old bineary-like better thickness - % cannot be expected - Ypp = min(Ypp,max(0.49,Ypps)); - Ypp = max(Ypp,min(0.51,Ypps)); - - % adaption only for very thin areas - Ypp(Ygmt<1.2) = Ypps(Ygmt<1.2); + % (*) EXTRA: magic sharpending function - enhancement of fine structures + if opt.sharpening + [Ypp, Ygmt] = sharpening(Ypp, Ygmt, vx_vol); end end % ====================================================================== -function [Ycd, Ywd] = cat_vol_cwdist(Yp0,opt) - - % CSF and WM distance - Ycd = zeros(size(Yp0),'single'); - Ywd = zeros(size(Yp0),'single'); +function Yp0 = WMTC(Yp0) +% WMTC. WM topology correction. - % multi-level distance estimation - hss = opt.levels; % number of opt.levels (as pairs) - for si = 1:hss - offset = max(0,min(.5, si * (.5)/hss)) / 2; - - Ycdl = cat_vbdist(single(Yp0 < ( 1.5 - offset)), Yp0 < 2.5 + 0.45*opt.extendedrange ); Ycdl(Ycdl > 1000) = 0; - Ycdh = cat_vbdist(single(Yp0 < ( 1.5 + offset)), Yp0 < 2.5 + 0.45*opt.extendedrange ); Ycdh(Ycdh > 1000) = 0; - - if opt.extendedrange - % additiona correction for levels>8 - Ycdlc = max(0,cat_vbdist(single(Yp0 < ( 2.5 - offset)), Ycdl>0 ) - .5 ); Ycdlc(Ycdlc > 1000) = 0; - Ycdhc = max(0,cat_vbdist(single(Yp0 < ( 2.5 + offset)), Ycdh>0 ) - .5 ); Ycdhc(Ycdhc > 1000) = 0; - Ycdl = Ycdl - Ycdlc; - Ycdh = Ycdh - Ycdhc; - end + numvx = 8; % size of correction + + % RD20231224: for WM topology smoothing + % the idea is to apply the voxel-based correction and _close_ small WM holes + % (<10 voxel) for the 2.75 and 2.25 boundaries for values above 2. - if opt.correctoffeset - if opt.correctoffeset==2 - offsetc = cat_stat_nanmedian(Ycdl(Ycdl>0) - Ycdh(Ycdl>0))/2; - else - offsetc = offset; - end - Ycdl(Ycdl>0) = max(eps, Ycdl(Ycdl>0) - (.5 - offsetc )); - Ycdh(Ycdh>0) = max(eps, Ycdh(Ycdh>0) + (.5 - offsetc )); - end + % indicate isolated holes and replace by median of the neighbors + Yp0(Yp0<0.35 & ~cat_vol_morph(Yp0<1,'l')) = 1; % close major wholes in the WM + Ymsk = Yp0==-1 & cat_vol_morph(Yp0>0.9,'d',1); % filter small wholes close to the WM + Yp0 = cat_vol_median3(single(Yp0)+1,Ymsk,~Ymsk)-1; - % mixing - Ycd = Ycd + .5/hss .* Ycdl + .5/hss .* Ycdh; - Ycw = max(0.1,.5/hss - opt.refine * Ycd/5 ); - %% - clear Ycdl Ycdh; + % indicate isolated objects and replace by median of the neighbors + Yp0(Yp0>0.65 & cat_vol_morph(Yp0==-1,'l')) = -1; + Ymsk = Yp0>0.95 & cat_vol_morph(Yp0<-0.9,'d',1); + Yp0 = cat_vol_median3(single(Yp0)+1,Ymsk,~Ymsk)-1; + % RD20210401: but there are some new background dots + Ymsk = Yp0==-1 & cat_vol_morph(Yp0>-1,'ldc',1.5); % close major wholes in the WM + Yp0 = cat_vol_median3(single(Yp0)+1,Ymsk,~Ymsk)-1; - % WM distances - Ywdl = cat_vbdist(single(Yp0 > ( 2.5 - offset)), Yp0 > 1.5 - 0.45*opt.extendedrange ); Ywdl(Ywdl > 1000) = 0; - Ywdh = cat_vbdist(single(Yp0 > ( 2.5 + offset)), Yp0 > 1.5 - 0.45*opt.extendedrange ); Ywdh(Ywdh > 1000) = 0; - if opt.extendedrange - Ywdlc = max(0,cat_vbdist(single(Yp0 > ( 1.5 - offset)), Yp0 > 1.05 ) - .5); Ywdlc(Ywdlc > 1000) = 0; - Ywdhc = max(0,cat_vbdist(single(Yp0 > ( 1.5 + offset)), Yp0 > 1.05 ) - .5); Ywdhc(Ywdhc > 1000) = 0; - Ywdl = Ywdl - Ywdlc; - Ywdh = Ywdh - Ywdhc; - end - if opt.correctoffeset - if opt.correctoffeset==2 - offsetc = median(Ywdl(Ywdl>0) - Ywdh(Ywdl>0))/2; - else - offsetc = offset; - end - Ywdl(Ywdl>0) = max(eps, Ywdl(Ywdl>0) + (.5 - offsetc )); - Ywdh(Ywdh>0) = max(eps, Ywdh(Ywdh>0) - (.5 - offsetc )); + % topology correction for WM based on the volume output of cat_vol_genus0 + % Closing of the WM object to avoid smaller holes. + % Opening is here not useful as the WM ! + tclevels = [2.75, 2.25]; + for ii = 1:1 % iterative ? + for tci = 1:numel(tclevels) + evalc(sprintf('Yppc = cat_vol_genus0(Yp0, %0.2f,0);',tclevels(tci))); + if 1 %tclevels(tci) > 2.5 % only closing + Yholes = Yppc==1 & Yp0=tclevels(tci); + Ybridges = Ybridges & ~cat_vol_morph(Ybridges,'l',[1e4,numvx]); % allow only small corrections + Yp0(Ybridges) = tclevels(tci) - 0.01; % open + clear Ybridges; end - - % mixing - Ywd = Ywd + Ycw .* Ywdl + (.5/hss*2 - Ycw) .* Ywdh; end + end + % object correction in the background (opening of GM-WM objects) + tclevels = [1.75, 1.25, 1.01]; + for tci = 1:numel(tclevels) + Yp0(Yp0>tclevels(tci) & ~cat_vol_morph(Yp0>tclevels(tci), ... + 'ldo',2.5 - tclevels(tci))) = tclevels(tci) - 0.01; + end end % ====================================================================== -function Yd = cat_vol_eudist(Yb,Ymsk,levels,correctoffeset) -%cat_vol_eudist. Euclidean distance estimation to mixed boundaries. -% ... +function [Ywd, Ycd, Ygmt] = keepdetails(Yp0, Ywd, Ycd, Ygmt, vx_vol, extendedrange,level) +% Although distances and thickness are quite good, PBT slightly trend to +% over-estimate the thickness in sulcal regions without CSF as the middle +% voxel is counted for both sides (simplified). Moreover, initial surface +% are partially created just on the original internal resolution (~1 mm), +% what result in blurring of small sucli. To keep these structures, we +% optimize regions where blurring/closing is happening by making them a +% bit thinner and correcting the CSF distance (keepdetails>0). The dual +% operation (of avoiding of blurring small gyri) can also be applied but +% is expected to have less effects as these structures are already quite +% save by the WM (keepdetails>1). % -% Yd = cat_vol_eudist(Yb,Ymsk,levels,correctoffeset) +% 1 -1.0*ppth ... 0 2.5235 ± 0.5959 mm 0.0686 / 0.0605 0.74% (9.16 mm²) 6.0 / 1.0 / 0.42% **** +% 3 -0.8*ppth ... 0 2.5403 ± 0.5872 mm 0.0682 / 0.0608 0.81% (10.01 mm²) 6.0 / 1.0 / 0.42% **** +% 2 -0.5*ppth ... 0 2.5671 ± 0.5727 mm 0.0677 / 0.0621 0.85% (10.49 mm²) 10.0 / 2.0 / 0.46% +% 2 -0.5*ppth ... .1 2.5515 ± 0.5863 mm 0.0690 / 0.0620 0.80% (9.87 mm²) 6.0 / 1.0 / 0.42% + + Ygmto = Ygmt; + + % estimate current percentage map (same as bellow) to identify + % and correct problematic areas + YM = Yp0 > 1.5 - extendedrange & ... + Yp0 < 2.5 + extendedrange & ... + Ygmt > eps; + Ycdc = Ycd; Ycdc(YM) = min(Ycd(YM), Ygmt(YM) - Ywd(YM)); + Ypp = zeros(size(Yp0),'single'); Ypp(Yp0 >= 2.5 + extendedrange) = 1; + Ypp(YM) = Ycdc(YM) ./ (Ygmt(YM) + eps); + % reduce thickness and CSF distance if this results in smoothing + % do this only in thin regions + for ppth = 1:-.1:0.1 + Yblurredsulci = Ypp=ppth,'c',1) & Ygmt < median(Ygmt(:))/2; + Yblurredsulci = cat_vol_smooth3X(Yblurredsulci, 2); + Ygmt = max( ... + max(0,Ygmto - 1/mean(vx_vol)), ... + Ygmt .* max(0.5,1 - 0.8 .* ppth*Yblurredsulci) ); + Ycd(Yblurredsulci>0) = max(0,min(Ycd(Yblurredsulci>0), Ygmt(Yblurredsulci>0) - Ywd(Yblurredsulci>0))); + end + + % reduce thickness and WM distance if this results in smoothing + % ... slow, worse values, + % - helpful to remove blue outliers? + if level == 2 + for ppth = 0.9:-0.1:0.1 + Yblurredgyri = Ypp>=ppth & cat_vol_morph(Ypp median(Ygmt(:))/2; + Yblurredgyri = cat_vol_smooth3X(Yblurredgyri, 2); + Ygmt = max( ... + min(0,Ygmto + 1/mean(vx_vol)), ... + Ygmt .* min(1.5,1 + .1 .* ppth*Yblurredgyri)); + Ywd(Yblurredgyri>0) = max(0,min(Ywd(Yblurredgyri>0), Ygmt(Yblurredgyri>0) + Ycd(Yblurredgyri>0))); + end + end +end +% ====================================================================== +function [Ypp, Ygmt] = sharpening(Ypp, Ygmt, vx_vol) +% sharpening. Reduce +% +% [Ypp, Ygmt] = sharpening(Ypp, Ygmt, vx_vol) +% -% add voxel size? + Ypp0 = Ypp; + + smoothfactor = .02; + smoothoffset = 1; + modthickness = 2; + + % create a sharpend version of the Ypps that enphalize regions that were smoothed + Ypps = Ypp + smoothfactor*1*(Ypp - cat_vol_smooth3X(Ypp,1/mean(vx_vol)) + smoothoffset * 0.000) + ... + smoothfactor*2*(Ypp - cat_vol_smooth3X(Ypp,2/mean(vx_vol)) + smoothoffset * 0.002) + ... + smoothfactor*4*(Ypp - cat_vol_smooth3X(Ypp,4/mean(vx_vol)) + smoothoffset * 0.004); + + % final smoothing to prepare surface reconstruction that also correct for WM topology issues + spm_smooth(Ypps,Ypps,.5/mean(vx_vol)); + Ypps = min(Ypp>0, max(0,Ypps)); % remove background artifacts + + % New version that only change mid-position values relevant for surface + % creation. Although this avoid the old binary-like better thickness + % cannot be expected + Ypp = min(Ypp,max(0.49,Ypps)); + Ypp = max(Ypp,min(0.51,Ypps)); + + % adopt thicnkess + Ymt = smooth3(Ypp - Ypp0) * modthickness; + Ygmt = Ygmt + Ymt; +end +% ====================================================================== +function [Yp0, Ywd, Ycd] = cleanupVessels(Yp0, Ywd, Ycd, doit) +% Filtering the Ywd to avoid mapping of artifacts by blood vessels or +% other uncleared tissues and regions. + + if doit + % == (1) cleanup distance values == + % estimate histogram + [hh,hr] = hist(Ywd(Ywd>eps) ./ Yp0(Ywd>eps),0:.1:50); + maxdist = max( hr( cumsum(hh)/sum(hh) < .99 ) ); + Ymsk = (Ywd ./ Yp0) > maxdist; + % correction + Yp0(Ymsk) = 1; + Ycd(Ymsk) = 0; + Ywd(Ymsk) = 0; + + + % == (2) cleanup by thickness values == + % estimate thickness + Ygmt = cat_vol_pbtp( round(Yp0) .* (Yp0 > 1.5) , Ywd, Ycd); + % remove outliers + gmtmn = cat_stat_nanmean(Ygmt(:)); + gmtsd = cat_stat_nanstd(Ygmt(:)); + Ymsk = Ywd < (Ygmt + gmtsd) & Ywd < (gmtmn + gmtsd); + Ygmt = Ygmt .* Ymsk; + gmtsd = cat_stat_nanstd(Ygmt(:)); + % approximate and strong smoothing + Ygmt = cat_vol_approx(Ygmt,'nh'); + Ygmt = cat_vol_smooth3X(Ygmt,16); + + % final correction + Ymsk = Ywd > (Ygmt + gmtsd) & Ywd > (gmtmn + gmtsd); + Ywd(Ymsk) = min(Ywd(Ymsk), max( Ygmt(Ymsk) + gmtsd , Ygmt(Ymsk)*1.25 )); + Yp0(Ymsk) = 1.5; + Ycd(Ymsk) = 0.5; + end +end +% ====================================================================== +function Yd = cat_vol_eudist(Yb, Ymsk, levels, correctoffeset) +%cat_vol_eudist. Euclidean distance estimation to mixed boundaries. +% +% Yd = cat_vol_eudist(Yb, Ymsk, levels, correctoffeset) +% +% Yd .. distance map +% Yb .. boundary map (with PVE) +% Ymsk .. masked regions for distance estimation +% +% levels .. number of dual distance measurements +% optimimum between 2 and 4 +% correctoffeset .. use generalized correction for the additional distance +% estimations, eg., for a more WM like value of 2.75 all +% distance values are assumed to be over +% (0 - none, 1 - default difference, 2 - estimated difference) +% +% see also cat_vol_pbtsimple:cat_vol_cwdist +% +% +% Todo: +% * Add voxel size? +% * Use as esternal function? +% if ~exist('Ymsk','var'), Ymsk = ~Yb; else, Ymsk = Ymsk>.5; end if ~exist('hss','var') @@ -438,3 +621,105 @@ Yd(Yd > max(size(Yb))) = 0; end +% ====================================================================== +function [Ycd, Ywd] = cat_vol_cwdist(Yp0,vx_vol,opt) +%cat_vol_cwdist. Estimation of CSF and WM distance in a label map Yp0. +% +% [Ycd, Ywd] = cat_vol_cwdist(Yp0,opt) +% +% Ycd, Ywd .. CSF and WM distance maps +% opt .. parameter structure +% .levels .. number of dual distance measurements +% .extendedrange .. estimate values also beyond the boundary to improve +% thickness mapping +% .correctoffeset .. use generalized correction for the additional distance +% estimations, eg., for a more WM like value of 2.75 all +% distance values are assumed to be over +% (0 - none, 1 - default difference, 2 - estimated difference) +% .range .. limitation to avoid bias by interpolation overshoot +% + + % CSF and WM distance + Ycd = zeros(size(Yp0),'single'); + Ywd = zeros(size(Yp0),'single'); + + opt.extendedrange = opt.extendedrange * opt.range; + + % additional correction map for values behind tissue boundary, e.g., + % for the WMD we stimate the distance from the GM/CSF boundary to + % limit WMD values to the maximal thickness value + % same idea as below + if opt.extendedrange + Ycdlc = Ycd; Ycdhc = Ycd; + Ywdlc = Ywd; Ywdhc = Ywd; + hss = opt.levels; % number of opt.levels (as pairs) + for si = 1:hss + offset = max(0,min(opt.range, opt.range * si/(hss+1))); + + Ycdlc = Ycdlc + 1/hss * max(0,cat_vbdist(single(Yp0 < ( 2.5 - offset)), Yp0 < 2.5 + opt.extendedrange ) - .5); + Ycdhc = Ycdhc + 1/hss * max(0,cat_vbdist(single(Yp0 < ( 2.5 + offset)), Yp0 < 2.5 + opt.extendedrange ) - .5); + + Ywdlc = Ywdlc + 1/hss * max(0,cat_vbdist(single(Yp0 > ( 1.5 - offset)), Yp0 > 1.5 - opt.extendedrange ) - .5); + Ywdhc = Ywdhc + 1/hss * max(0,cat_vbdist(single(Yp0 > ( 1.5 + offset)), Yp0 > 1.5 - opt.extendedrange ) - .5); + end + Ycdlc(Ycdlc > 1000) = 0; Ycdhc(Ycdhc > 1000) = 0; + Ywdlc(Ywdlc > 1000) = 0; Ywdhc(Ywdhc > 1000) = 0; + end + + + % multi-level distance estimation + hss = opt.levels; % number of opt.levels (as pairs) + for si = 1:hss + offset = max(0,min(opt.range, opt.range * si/(hss+1))); + + Ycdl = cat_vbdist(single(Yp0 < ( 1.5 - offset)), Yp0 < 2.5 + opt.extendedrange ); Ycdl(Ycdl > 1000) = 0; + Ycdh = cat_vbdist(single(Yp0 < ( 1.5 + offset)), Yp0 < 2.5 + opt.extendedrange ); Ycdh(Ycdh > 1000) = 0; + + if opt.extendedrange + Ycdl = Ycdl - Ycdlc; + Ycdh = Ycdh - Ycdhc; + end + + if opt.extendedrange + if opt.correctoffeset==2 + offsetc = offset/mean(vx_vol) + (cat_stat_nanmedian(Ycdl(Ycdl>0 & Ycdh>0) - Ycdh(Ycdl>0 & Ycdh>0)))/2; + else + offsetc = offset/mean(vx_vol); + end + Ycdl(Ycdl>0) = max(eps, Ycdl(Ycdl>0) - (.5 + offsetc) ); + Ycdh(Ycdh>0) = max(eps, Ycdh(Ycdh>0) + (.5 + offsetc) ); + end + + % mixing + Ycd = Ycd + .5/hss .* Ycdl + .5/hss .* Ycdh; + %Ycw = max(0.1,.5/hss - opt.refine * Ycd/5 ); + % idea was to could the boundaries different depending on the CSF distance + clear Ycdl Ycdh; + + + % WM distances + Ywdl = cat_vbdist(single(Yp0 > ( 2.5 - offset)), Yp0 > 1.5 - opt.extendedrange ); Ywdl(Ywdl > 1000) = 0; + Ywdh = cat_vbdist(single(Yp0 > ( 2.5 + offset)), Yp0 > 1.5 - opt.extendedrange ); Ywdh(Ywdh > 1000) = 0; + + if opt.extendedrange + Ywdl = Ywdl - Ywdlc; + Ywdh = Ywdh - Ywdhc; + end + if opt.correctoffeset + if opt.correctoffeset==2 + offsetc = offset + (cat_stat_nanmedian(Ywdl(Ywdl>0 & Ywdh>0) - Ywdh(Ywdl>0 & Ywdh>0)))/2; + else + offsetc = offset; + end + Ywdl(Ywdl>0) = max(eps, Ywdl(Ywdl>0) + (.5 + offsetc) ); + Ywdh(Ywdh>0) = max(eps, Ywdh(Ywdh>0) - (.5 + offsetc) ); + end + + % mixing + Ywd = Ywd + .5/hss .* Ywdl + .5/hss .* Ywdh; + %Ywd = Ywd + Ycw .* Ywdl + (1/hss - Ycw) .* Ywdh; + end + + +end +% ====================================================================== diff --git a/cat_vol_set_com.m b/cat_vol_set_com.m index d7fb60f9..f115b51a 100644 --- a/cat_vol_set_com.m +++ b/cat_vol_set_com.m @@ -36,7 +36,12 @@ % pre-estimated COM of MNI template com_reference = [0 -20 -15]; -fprintf('Correct center-of-mass \n'); +if nargin == 1 + % call from cat_run_job that will add the time and then the line break + fprintf('Correct center-of-mass '); +else + fprintf('Correct center-of-mass \n'); +end for i=1:n Affine = eye(4); if isfield(V(i),'dat')