1) New: Added final CSV report to cat_run with changes in cat_io_xml2…

…csv.

			   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].
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

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:

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','.'))
 

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

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

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

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<ASGLU>
           catch apperr
             %% very simple affine preprocessing ... only simple warning

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<ASGLU>
           catch apperr
             %% very simple affine preprocessing ... only simple warning

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 && SIs<SIOs*0.9 && iter<maxiter
+      SIOs = 100; SIs = 80; maxiter = 1; iter = 0; 
+      while SIs>5 && SIs<SIOs*0.9 && iter<=maxiter
         SIOs = SIs; iter = iter + 1; 
-        [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)); 
+        [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');