From 102cca711de44510cee6ee549607b5a131258f5d Mon Sep 17 00:00:00 2001
From: Christian Gaser <christian.gaser@uni-jena.de>
Date: Tue, 21 May 2024 01:23:33 +0100
Subject: [PATCH] Fixed: Corrected bug in cat_vol_slice_overlay, where clipping
 with OV.func was        not correct and atlas regions were wrong in that
 case.

Changed paths:
 M CHANGES.txt
 M cat_main_LAS.m
 M cat_vol_iscale.m
 M cat_vol_slice_overlay.m
---
 CHANGES.txt             | 11 +++++++++++
 cat_main_LAS.m          |  2 +-
 cat_vol_iscale.m        |  4 ++--
 cat_vol_slice_overlay.m | 30 +++++++++++++++++-------------
 4 files changed, 31 insertions(+), 16 deletions(-)

diff --git a/CHANGES.txt b/CHANGES.txt
index 5dc37abf..43e993c1 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -1,4 +1,15 @@
 ------------------------------------------------------------------------
+r2582 | gaser | 2024-05-21 01:23:31
+
+Changed paths:
+ M CHANGES.txt
+ M cat_main_LAS.m
+ M cat_vol_iscale.m
+ M cat_vol_slice_overlay.m
+
+Fixed: Corrected bug in cat_vol_slice_overlay, where clipping with OV.func was 
+       not correct and atlas regions were wrong in that case.
+------------------------------------------------------------------------
 r2579 | gaser | 2024-05-10 21:01:35
 
 Changed paths:
diff --git a/cat_main_LAS.m b/cat_main_LAS.m
index cd63e8af..cf3ee81c 100644
--- a/cat_main_LAS.m
+++ b/cat_main_LAS.m
@@ -323,7 +323,7 @@
   %            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
+  %    Ysw4 .. was removed 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
diff --git a/cat_vol_iscale.m b/cat_vol_iscale.m
index da598327..dfd9adc7 100644
--- a/cat_vol_iscale.m
+++ b/cat_vol_iscale.m
@@ -1,5 +1,5 @@
 function [TI,varargout] = cat_vol_iscale(T,action,vx_vol,varargin)
-% CAT Prprocessing Intensity Scaling Functions
+% CAT Preprocessing Intensity Scaling Functions
 % ______________________________________________________________________
 % Set of functions for intensity scaling of an image T. 
 %
@@ -213,7 +213,7 @@
         WIr  = cat_vol_resize(WIrr,'dereduceV',resT2); clear WIrr;
         WIr  = WIr * median(TSr(M2r(:) & Gr(:)<Gth & TSr(:)>WIr(:)*0.95) ./ WIr(M2r(:) & Gr(:)<Gth & TSr(:)>WIr(:)*0.95));
 
-        %% rough skull-stipping
+        %% rough skull-stripping
         % use the head distance map D to find a central tissue (low gradient) with WM intensity
         [Grr,Trr,GWMrr,Drr,M2rr,resT4] = cat_vol_resize({Gr,TSr./WIr,GWMr./WIr,Dr,M2r},'reduceV',vx_vol,2,32,'meanm');
         Brr  = cat_vol_morph(Trr<1.2 & Trr>max(0.8,GWMrr) & Grr<Gth & M2rr); 
diff --git a/cat_vol_slice_overlay.m b/cat_vol_slice_overlay.m
index fc23a6e6..a16073bd 100644
--- a/cat_vol_slice_overlay.m
+++ b/cat_vol_slice_overlay.m
@@ -185,11 +185,12 @@
     slices{1} = OV.slices;
   else
     SO.img(2).vol = spm_vol(OV.name);
-    [mx, mn, XYZ, img2] = volmaxmin(SO.img(2).vol);
+    if isfield(OV,'func')
+      [mx, mn, XYZ, img2] = volmaxmin(SO.img(2).vol,OV.func);
+    else
+      [mx, mn, XYZ, img2] = volmaxmin(SO.img(2).vol);
+    end
 
-    % apply function to img2
-    if isfield(OV,'func'), i1 = img2; eval(OV.func); img2 = i1; end
-    
     % threshold map and restrict coordinates
     Q = find(compare_to_threshold(img2,range(1)) & le(img2,range(2)));
     XYZ = XYZ(:, Q);
@@ -259,7 +260,11 @@
 SO.slices = slices;
 
 if isempty(SO.slices)
-  [mx, mn, XYZ, vol] = volmaxmin(SO.img(2).vol);
+  if isfield(OV,'func')
+    [mx, mn, XYZ, vol] = volmaxmin(SO.img(2).vol,OV.func);
+  else
+    [mx, mn, XYZ, vol] = volmaxmin(SO.img(2).vol);
+  end
   
   % threshold map and restrict coordinates
   Q = find(compare_to_threshold(vol,SO.img(2).range(1)));
@@ -544,7 +549,7 @@
 
 % atlas labeling
 if ~isempty(xA)
-  [mx, mn, XYZ, vol] = volmaxmin(SO.img(2).vol);
+  [mx, mn, XYZ, vol] = volmaxmin(SO.img(2).vol,SO.img(2).func);
   
   % threshold map and restrict coordinates
   if SO.img(2).range(1) >= 0
@@ -555,9 +560,7 @@
   M = SO.img(2).vol.mat;
   XYZmm = M(1:3, :) * [XYZ; ones(1, size(XYZ, 2))];
   
-  % apply func that is defined for "i1"
   i1 = vol;
-  eval(SO.img(2).func)
   
   % remove NaN values
   Q = find(isfinite(i1));
@@ -794,7 +797,7 @@
 return
 
 % --------------------------------------------------------------------------
-function [mx, mn, XYZ, img] = volmaxmin(vol)
+function [mx, mn, XYZ, img] = volmaxmin(vol,func)
 
 if nargout > 2
   XYZ = [];
@@ -805,11 +808,12 @@
 
 mx = -Inf; mn = Inf;
 for i = 1:vol.dim(3)
-  tmp = spm_slice_vol(vol, spm_matrix([0 0 i]), vol.dim(1:2), [0 NaN]);
-  tmp1 = tmp(isfinite(tmp(:)) & (tmp(:) ~= 0));
+  i1 = spm_slice_vol(vol, spm_matrix([0 0 i]), vol.dim(1:2), [0 NaN]);
+  if nargin > 1, eval(func); end
+  tmp1 = i1(isfinite(i1(:)) & (i1(:) ~= 0));
   if ~isempty(tmp1)
     if nargout > 2
-      [Qc Qr] = find(isfinite(tmp) & (tmp ~= 0));
+      [Qc Qr] = find(isfinite(i1) & (i1 ~= 0));
       if size(Qc, 1)
         XYZ = [XYZ; [Qc Qr i * ones(size(Qc))]];
         if nargout > 3
@@ -879,7 +883,7 @@
   if i == 1
     SO.img(i).cmap = gray;
     %[mx, mn] = volmaxmin(SO.img(i).vol);
-    [tmp, th]=cat_stat_histth(spm_read_vols(SO.img(i).vol),0.95,0);
+    [tmp, th] = cat_stat_histth(spm_read_vols(SO.img(i).vol),0.95,0);
 
     SO.img(i).range = th;
   else